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1 .   INTRODUCTION 

Define  an  automatic  program  as.  one  that  can 

1.  take  as  input  members  of  some  broad  class  of  problems, 

2.  try  to  control  errors  within  specified  limits, 

3.  try  to  minimize  the  work  required  for  the  solution. 

The  problem  of  the  present  paper  is  to  look  for  an  automatic 
method  that  can  be  applied  to  some  subset  of  partial  differential  equations. 
The  primary  value  of  a  program  involving  an  automatic  method  would  be  to 
provide  the  person  who  needs  only  a  few  runs  with  a  particular  type  of 
equation  the  opportunity  to  find  a  solution,  within  acceptable  tolerances, 
without  having  to  write  and  debug  a  specialized  program.   In  this  environ- 
ment, the  program  will  have  no  information  about  the  solution,  or  how  errors 
might  accumulate.  While  maintaining  sufficiently  small  errors,  an  attempt 
needs  to  be  made  to  keep  the  cost  of  such  a  solution  as  low  as  possible.   To 
accomplish  these  goals,  decisions  need  to  be  made  during  the  integration  for 
balancing  the  two  conditions  above  based  on  some  sort  of  error  estimate  and 
the  user's  error  tolerance. 

In  the  present  study,  the  subclass  of  partial  differential  equations 
was  chosen  to  be  the  set  of  second  order  mixed  initial  boundary  value  problems 
in  two  variables  (x  and  t)  of  the  form  F(u,  u  ,  u   ,  x,  t)  =  u  .  The  general 

X     XX  T> 

method  used  involves  integrating  along  lines  which  are  parallel  to  the  time 
axis  and  which  pass  through  a  finite  set  of  points  on  the  spatial  axis.  To 
obtain  an  idea  of  this  method,  consider  the  equation 

9u  _  afu 
9t  ~   2 

dX 


with  initial  condition 

u(x,0)  =  sin(Trx)  x  in  [0,l] 
and  boundary  conditions 

u(0,t)  =u(l,t)  =  0 
To  begin,  select  a  partition 

P  =  {0  =  xQ,  x^  x2,..  .,  xn_1,  xn  =  1} 

and  then  define  a  set  of  functions  u. (t)  such  that 

u.(0)  =  u(x. ,0)   i  =  1,  2,...,  n. 

By  the  use  of  some  interpolating  polynomial,  a  discrete  analog  is  obtained 
for  the  spatial  derivative,  which  will  be  represented  by 


in 

8x2 


X. 

1 


2  W 


Using  this  analog,  we  obtain  the  discrete  system 


du. 

-~  =  A  (u.  )   i  =  1,.  ..  ,  n 
dt     n  i 


This  system  differs  from  the  system 


9u(x,t) 


3t 


3  u(x,t) 


3x 


x. 


which  is  obtained  from  the  true  solutions  by  an  amount  6(i,t). 

The  analog  A  (u. )  and  the  partition  P  are  chosen  in  an  attempt 
n  i 

to  keep  u. (t)  -  u(x. ,  t)  within  the  user's  desired  tolerance. 

Once  the  partition  and  analog  are  chosen,  one  then  integrates  the 
system  in  the  time  direction  along  the  lines  implied  by  P  from  t  =  t  to 


t  =  t 


If  at  any  time  during  this  integration  it  appears  that  the  error 


is  too  large,  the  step  is  repeated  with  the  possibility  of  changing  either 
AQ(ui)  or  P.  A  change  in  the  partition  would  mean  continuing  the  integra- 
tion along  a  new  set  of  lines. 

There  is  also  the  possibility  that  the  error  is  much  smaller  than 
is  required.  Again,  modifications  in  A^iu)  or  P  might  be  attempted  so 
that  the  work  of  obtaining  the  desired  error  is  minimized.   These  decisions 
will  need  to  be  made  on  the  basis  of  some  estimate  of  the  error. 

The  basis  for  the  integration  of  the  ordinary  differential  equation 
system  above  is  a  program  (DIFSUB)  developed  by  Gear  [l]  for  automatically 
solving  systems  of  stiff  ordinary  differential  equations.  As  an  automatic 
program,  it  attempts  to  keep  the  errors  introduced  in  integrating  in  the  t 
direction  within  a  given  tolerance . 

In  integrating  along  the  possibly  changing  set  of  lines,  there  are 
thus  three  ways  of  trying  to  keep  the  errors  in  the  solution  at  a  desired 
value.   First,  the  partition  can  be  changed.  Second,  the  nature  of  the  analog 
used  to  discretize  the  system  can  be  altered.   And  third,  the  amount  of  the 
error  parameter  given  to  the  integration  routine  can  be  modified. 

In  order  to  achieve  an  automatic  program  for  the  entire  system, 
several  items  will  need  to  be  considered: 

1.  the  forms  of  the  analogs  to  be  used, 

2.  the  methods  of  obtaining  error  estimates, 

3.  the  methods  of  obtaining  and  modifying  partitions, 

h.      the  implementation  of  the  user's  equations  in  an  easily 
usable  form. 


2 .   PROCEDURES 

2.1  Spatial  Discretization 

There  are  several  problems  to  be  considered  in  looking  for  a 
spatial  analog.   These  include  the  number  of  equations  required,  the  ease 
of  changing  the  partition,  and  the  situation  at  the  boundaries.   Three 
methods  were  initially  considered: 

1.  Lagrangian  methods, 

2.  Hermite  type  methods, 

3.  Linear  transformations  of  the  Lagrangian  methods. 
These  will  then  be  used  on  equations  of  the  form 

—  =  f(u,  u  ,  u   ,  x,  t) 
3t         x   xx 

The  equations  are  then  put  into  a  discrete  form  by  replacing  u  by  a  set 

of  functions  u. (t)  and  replacing  the  right  side  by  an  analog  A  (u. ) .  This 
l  n  i 

then  gives  the  system  of  equations 

du. 

ar=VV  w 

In  the  case  of  the  Lagrangian  polynomials,  there  are  several 

possibilities.   First,  there  is  the  class  of  polynomials  using  equidistant 

points.   If  the  number  of  points  (n)  in  the  interpolating  polynomial  is  odd, 

then  it  is  possible  to  take  advantage  of  symmetry  properties  to  eliminate 

one  of  the  terms  in  the  error.   Given  the  approximation 

n 
y(x)  =  £  I    (x)  f(a  )  (2) 

i=l 

to  a  function  f(x)  for  interpolating  points  a.,...,  a  ,  the  error  term  is 

given  by 


-Sr-  *m 

where  p^x)  =  (x-^ ) . . .  ( x-aQ )  and  n  is  some  point  in  the  interval  containing 
these  points.   Taking  the  derivative  of  the  interpolating  polynomial,  an 
approximation  is  obtained  for  f'(.x). 


n 


y'Cx)  =  z  i  ;Cx)  f(a.)  to) 

1=1  x  x 

To   obtain  ideas   about  the  error,   take  the  derivative  of  the  error  term  of 
(2)    (see  2.2   for  discussion  and  development).      For  equation    (3),   the  error 
term  would  be  given  by 

n!        f       (ril)+InTnTf  V  (4) 

Since  pn(x)   =  0  at   tabular  points,  where   all  such  evaluations  will  be  made, 
(k)   reduces  to 

_n An),       v 

^T^f       (rii}  (5) 

If  this  process  is  repeated,  the  error  for  the  second  derivative  is 

^f(n»K)+^f(n+1)(n2)  +  ^f(n+2)(n3)  (g) 

If  the  interpolating  polynomial  uses   an  odd  number  of  points,    and  if  the 
derivative   is  estimated  for  the  center  point,    then  the   first   term  is    zero. 
Consider 

n        n              Pn(x> 
P"(x)    =     I        I      , "L 

i=l  j=l   (x-a^tx-a   ) 

Since  evaluation  is  at  the  center  point,  x  -  -^  =  0 ,  and  all  terms 
containing  this  factor  will  vanish.   As  a  result,  P;'(x)  can  be  rewritten  as 

(where  the  subscript  nl  is  2^) 

2  ' 


n  2-p  U) 

iftexl 

Further,  since  the  points  are  symmetrically  located  around  a 


CT) 


nl 


(a   -a   +i)=(an-a,-i).  Noting  this  symmetry,  (T)  can 
nl    nl         nl    nl 


he  written 


hi        2-p  (x)  2'pJx) 


PnU)  =  !=1  '^V^V^TT  +  U-anl)U-anl  =  D 

=  0 
Since  the  last  term  of  (6)  is  also  zero,  the  expression  reduces  to 

2p'(x)  ,    ._  v 

^TT^^In)  (8) 

Since  p  (x)  is  a  function  of  the  spacing  h,  the  order  could  he  defined 
as  the  lowest  power  of  h  contained  in  the  error.   In  (5)  and  (8)  these 
orders  are  the  same,  so  that  the  truncation  errors  should  react  similarly 
to  changes  of  h.   In  the  event  the  symmetry  is  lost,  the  order  associated 
with  the  second  derivative  is  reduced  hy  one,  and  P^(x)  is  no  longer  zero. 
The  error  term  for  the  second  derivative  would  prohably  then  he  the 
determining  factor  in  selection  of  an  h. 

One  of  the  main  difficulties  with  this  class  of  methods  results 
from  the  equidistant  spacing.   If  the  solution  is  "non -smooth" ,  the 
spacing  for  the  partition  will  have  to  he  set  so  as  to  reduce  the  error  at 
the  worst  point.  This  could  lead  to  a  much  larger  number  of  equations  than 
needed,  for  in  the  smoother  regions  the  error  would  he  expected  to  be  much 
less  than  is  required.   Ideally,  one  is  aiming  to  have  partition  points 


7 

such  that  the  total  error  introduced  at  each  is  at  the  level  desired  by 
the  user. 

In  an  attempt  to  keep  the  error  closer  to  the  user's  request, 
one  could  divide  the  interval  into  sub intervals  with  spacings  that  are 
locally  equidistant.   If  the  analog  uses  five  or  more  points,  such 
subdivisions  have  problems  near  the  end  of  each  subinterval.   The 
derivative  values  near  these  subinterval  end  points  would  either  have  to 
use  a  nonsymmetric  form  or  a  point  from  an  adjacent  interval  with  possibly 
a  different  step  size. 

When  using  a  three-point  analog,  these  problems  are  reduced. 
Each  analog  could  use  spacings  that  are  multiples  of  some  h.   In  order 
to  circumvent  the  problem  of  the  end  point,  one  could  consider  it  as  the 
center  point  of  another  subinterval  overlapping  those  on  either  side. 
For  example,  given  the  partition 


p6 


P3PUP5 


and  an  h  ' 


15  3    71 
0  ITIS8I6T 

rr,  the  analogs  could  be  set  up  in  the  following  way 
derivative  at  points  used 


vv  y2,   p6 

P2>  p3'  Vh 
P3»  P^,  P5 

vk,  p5,  p6 

Pr  P6,  P? 


8 

In  practice,  enough  points  would  be  added  so  that  several  of  the  sub- 
intervals  would  not  be  contained  within  another  interval. 

Another  way  to  try  to  keep  all  errors  somewhat  close  to  the 
user's  request  would  be  to  use  unequal  spacings  throughout  the  interval. 
There  are  several  disadvantages  to  this  approach.  First,  the  order  of 
the  error  term  would  be  lower  than  those  associated  with  symmetric- 
equispaced  ones.   Second,  there  exists  the  problem  of  calculating  the 
derivatives.  With  equidistant  points,  a  set  dependence  exists  between 
values  of  the  function  and  the  derivative  analog  at  any  partition  point . 
In  the  case  of  unequal  spacings,  either  the  analog  would  have  to  be 

calculated  each  time  it  was.  needed  or  the  dependence  of  A  (u.  )  on  each 

n  1 

point  would  have  to  be  calculated  and  then  stored  for  a  given  partition. 
Third,  the  initial  partition  would  be  somewhat  more  difficult  to  develop 
because  of  the  freedom  of  changing  the  spacings  between  the  points.  As 
a  result,  it  could  happen  that  the  spacings  between  any  two  pairs  of 
adjacent  points  could  be  greatly  different.   In  this  case,  it  would  be 
desirable  to  go  back  and  change  other  points  so  as  to  avoid  widely 
differing  spacings  for  neighboring  points. 

A  second  class  of  analogs  involves  those  methods  that  use  more 
than  just  the  functional  value  at  each  point.   The  primary  subclass  would 
be  the  Hermite  interpolating  polynomials.  With  polynomials  of  this  form, 
it  is  necessary  to  save  the  first  derivative  at  each  point  as  well  as  the 
functional  value.   In  order  to  develop  this  additional  value  at  each  time 
step,  it  would  be  necessary  to  solve  an  additional  equation  at  each  point. 
Such  an  equation  can  be  found  by  taking  the  partial  derivative  with  respect 
to  the  space  variable  of  both  sides  of  the  given  equation 


Ut  =  f(u'  V  Uxx'  x'  t)  (9) 

This  gives 

3lUt  =  3x"f(u'  V  uxx'  X'  t}  ("J 

Since  the  solution  u  is  assumed  to  be  continuous  and  to  have  as  many- 
continuous  partials  as  may  be  required,  the  order  of  differentiation  is 
changed  to  give 

3t  Ux  =  ^f(u'  V  Uxx>  x> t}  (id 


At  this  point,  u  is  replaced  by  a  new  variable  v.  and  the 
x  l 

right  side  of  (11)  is  replaced  by  a  discrete  analog  B  (u. ,  v.),  giving 


n'  i'   i 


3 


3?Vi  =  VV  vi}  (12) 

The  solution  v.  of  this  addition  equation  is  not  expected  to  be  ill-behaved 
provided  the  solution  is  continuous.   Consider  the  Taylor  series  expansion 
of  u 

u(x,t)  =  u(x0,tQ)  +  ux(v0,t0)(x-x0)  +  ut(x0,t0)(t-tQ)  +.,.. 

The  presence  of  ux  in  this  expansion  indicates  that  any  unusual  behavior 
that  might  be  reflected  in  ux  would  be  present  in  u.   As  a  result,  if  the 
original  system  of  equations  (9)  can  be  solved,  then  it  is  expected  that  the 
combined  system  of  (9)  and  (ll)  can  also  be  solved. 

Just  as  in  the  Lagrangian  case,  there  is  an  advantage  in  using 
the  symmetry  properties  to  eliminate  terms  in  the  expression  for  the  error. 
Further,  all  of  the  problems  associated  with  the  Lagrangian  polynomials — 
when  more  than  three  points  are  involved— are  compounded  by  the  fact  that 
the  first  derivative  must  be  calculated  and  saved  for  each  point.   A  new 
problem  exists  in  that  additional  values  are  needed  at  the  boundary.   In 


10 

the  case  of  the  original  equations,  the  boundary  conditions  provide 

additional  equations  that  can  he  solved  with  the  differential  equations . 

With  Vn  and  v  at  the  boundaries,  it  will  be  necessary  to  develop  other 

equations,  or  at  least  provide  values,  for  these  variables. 

Another  possible  member  of  this  second  class  of  interpolating 

polynomials  would  be  those  that  match  first  and  second  derivatives  as 

well  as  the  functional  values  at  the  interpolating  points .   The  use  of 

such  polynomials  (trimite)  and  the  problems  involved  will  be  very  similar 

to  those  encountered  in  the  Hermite  case.   Here,  two  additional  equations 

will  need  to  be  supplied  for  each,  point.   The  first  of  these  would  be  the 

same  as  the  additional  equation  in  the  Hermite  case.   The  second  would  be 

found  by  differentiating  both  sides  of  (ll)  with  respect  to  x  giving 

2 
3   9       3   „/  .  \ 

T"  TT  u    — T  f(u,  u  ,  x   ,  x  t) 
3x  3t  x   _  2       x   xx   , 
3x  ' 

dv. 

Interchanging  the  order  again,  introducing  an  additional  variable  w.  =— — 

l   dx 

and  an  analog  C  (u. ,  v. ,  w. ) ,  the  discrete  version  becomes 
n  l   l   l 

dw. 

jr=cn(v  vwi>  (13> 

Another  problem,  that  was  also  present  in  the  Hermite  but  becomes 
more  evident  in  the  trimite  polynomials,  is  the  presence  of  higher  ordered 
derivatives.   The  error  term  for  the  trimite  will  contain  the  ninth  and 
tenth  derivatives.   Although  accuracy  will  not  be  required,  a  reasonable 
approximation  is  desirable. 

The  one  additional  possibility  that  was  considered  for  obtaining 
an  analog  was  to  use  a  linear  transformation  of  a  Lagrangian  method  which  would 
save  different  quantities,  but  which  would  possibly  save  work.   Suppose 


11 


eight  interior  points  and  two  boundary  points  (0  =  x  ,...,  x  =  1}  are 
used  with  the  equation  ut  =  u^  x  e[0,l].   The  discretized  equations 


du 


would  be  represented  by  _£  =  Au  where  A  is  given  by 


Q 

0 

Q 

0 

0 

0 

0 

0 

0 

0 

1 

-2 

1 

0 

0 

0 

0 

0 

0 

0 

Q 

1 

-2 

1 

0 

0 

0 

0 

0 

0 

0 

0 

1 

-2 

1 

0 

0 

0 

0 

0 

0 

0 

0 

1 

-2 

1 

0 

0 

0 

0 

0 

0 

0 

0 

1 

-2 

1 

0 

0 

0 

Q 

0 

0 

0 

0 

1 

-2 

1 

0 

0 

0 

0 

0 

0 

0 

0 

1 

-2 

1 

0 

0 

0 

0 

0 

0 

0 

0 

1 

-2 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

n, 

it  would  be 

easier  t< 

3  save 

l   u_  . 

l]     .     l 

i        i 

l"    » 

1'     l'  "3  s  "3'  "5'  u5' 
uT,  u^  instead  of  the  functional  values  at  the  eight  interior  points. 

The  transformation  T  would  then  be  given  by 


1 

0 

0 

0 

0 

0 

0 

0 

0 

0 
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1 
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0 

0 

0 

1 

-2 

1 

0 

0 

0 

0 

0 

0 

0 

0 

1 

-2 

1 

0 

0 

0 

0 

0 

0 

0 

0 

1 

-2 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

1 

-2 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

1 

-2 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

12 


-1 


The  transformed  matrix  B  =  TAT    "  then  become 


0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

2 

-3 

0 

l 

a 

0 

0 

0 

0 

0 

0 

0 

0 

0 

l 

0 

0 

0 

0 

0 

a 

1 

0 

-2 

Q 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

1 

0 

-2 

0 

1 

0 

0 

0 

0 

a 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

1 

0 

-2 

0 

1 

0 

0 

0 
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0 

0 

0 

0 

0 

0 

J 


2 
For  3u  =  l_u  the  transformation  gives  results  that  would  perhaps  he 


3t 


3x 


worth  using.   However,  if  the  original  equation  contained  the  first 
derivative  as  well  as  the  second,  a  second  matrix  D  =  TCT-1  needs  to 
he  calculated  where  2'C  is 


0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

-1 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

-1 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

-1 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

-1 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

-1 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

-1 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

-1 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

-1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

L' 
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and  2*D  is 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

i    2 

-2 

-1 

0 

0 

0 

0 

0 

0 

0 

-k 

5 

2 

-1 

0 

0 

0 

0 

0 

0 

-2 

4 

2 

-2 

-1 

0 

0 

0 

0 

0 

k 

-7 

-1+ 

1+ 

2 

-1 

0 

0 

0 

0 

2 

-It 

_P 

k 

2 

-2 

-1 

0 

0 

0 

-1+ 

8 

It 

-1 

-It 

It 

2 

-1 

0 

0 

-2 

U 

2 

-k 

-2 

1+ 

2 

-2 

-1 

0 

1+ 

-8 

J* 

8 

1+ 

-7 

-1+ 

It 

2 

-1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

Since  the  user's  equation  is  unknown,  the  need  for  both  matrices  is  extremely 
likely.   Therefore,  the  idea  of  saving  the  second  derivative  loses  its  appeal, 

If  the  function  and  the  first  derivative  are  saved,  then  the 
results  are  similar.   For  the  first  derivative,  the  resulting  matrix  is 


(B  = 

TAT"1 ) . 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

1 

"It 

0 

1 
It 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

B  = 

0 

1 
IT 

0 

1 

2 

0 

1 
IT 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

1 
IT 

0 

1 

~2 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

1 

it 

0 

1 

2 

0 

1 
IT 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

lU 


,-!• 


and  that  for  the  second  derivative  is  (D  =  TCT   ) 


D  = 


0 

0 

0 

0 

0 

0 

G 

0 

a 

0 

2 

-2 

-2 

0 

Q 

0 

Q 

0 

0 

0 

1 

1 
2 

2 

1 

2 

0 

0 

0 

0 

0 

0 

2 

0 

-k 

-2 

-2 

0 

0 

0 

0 

0 

0 

1 
2 

0 

0 

1 

"2 

0 

0 

0 

0 

2 

0 

-k 

1 
2 

-u 

-2 

-2 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-2 

1 
2 

0 

0 

2 

0 

-i+ 

a 

-k 

1 
2 

-u 

-2 

-2 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-2 

1 

2 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

Since  the  problems  the  procedures  are  meant  to  solve  could  involve  both 
derivatives,  the  -whole  idea  of  using  linear  transformations  does  not 
appear  to  provide  anything  of  interest. 

As  a  result  of  the  disadvantages  of  going  to  more  than  three  points 
it  was  decided  that  three  methods  would  be  used.   These  were  the  Lagrangian, 
Hermite,  and  trimite  polynomials  using  three  points.   This  would  provide 
methods  of  three  different  orders.  Although  the  Hermite  and  trimite  forms 
have  additional  problems,  the  flexibility  of  having  the  different  ordered 
methods  and  the  possibility  of  a  fewer  number  of  equations  makes  consideratic 
of  these  methods  worth  undertaking. 

If  the  polynomials  for  Lagrangian,  Hermite,  and  trimite  are  called 


Y  ,  Y  ,  Y  ,  respectively,  the  resulting  formulas  can  be  written 
L   H   T 
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n 


Y     =     I     Z.(x)    f(x.) 
L       i=l     X  X 


n 


YH  = 


=     L      {h  (x)    f(x.  )   +  h,(x)    f'(x.) 

i-1  -L  1  1  1 


n 


T  -     2      (t    (x)    f(x   )   +  t   (x)    f»(x.)   +  t.(x)    f"(x.)) 
i=l  j.  x  ii  x 

where  t.(x),   £(x),   and  t±(x)   are  polynomials   associated  with  the  trimite 
method. 

The  development  of  t,   t,   t  parallels  the   development  of  the 
Hermite  polynomials. 
Let 

\  =    (aik   ^  +   a2k  X  +   V  *k(x) 

\  =    (blk  ^  +  a2k  x  +  a3k}  *k(x) 

\  =    (cik  ^  +  C2k  X  +  C3k}  £k(x) 
There  are  three  conditions  for  each  of  the  equations 

?k(Xj )  =  o    t.(Xj )  =  o    !»(«, >  .  «k. 

Solving  each  row  gives 
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blk  ■  -*£<v 

1,2k-1  +  6ai1i0,k) 

V=  "Xk(1  +  3xk£i(xk)) 

_  I 

°lk  "  2 

C2k  =  -\ 

1  2 
c3k  "  2*k 

Substituting  into  the  original  expressions  and  simplifying,  the  equations 
become 

tk(x)  =  [i  -  §<*-V2  *'i(V  +  6(x-\)2  l't(\))2  "  3<x-\'- 

tk(x)=  [(x-^)  -3tk(xk)(x-xk)2l^(x) 
tk(x)  -  §<x-x/  f3(x) 

Using  a  proof  by  Wendroff  [2],  the  error  term  is 

E^-f(n)(n)  . 
n! 

where  n  is  the  number  of  values  involved.   Since  a  three-point  formula  for  the 
derivative  of  the  center  point  is  being  used,  the  coefficients  for  each  of 
these  interpolating  polynomials  can  be  determined. 


IT 

When  using  the  Lagrangian  polynomials ,  the  given  equations  are 
used,  so  the  first  and  second  derivatives  are  needed.   These  are 

YjJ  =  (f(x+h)  -  f(x-h))/2h 

Y£   =  (f(x+h)  -  2f(x)  +  f(x-h))/h2 

In  the  case  of  the  Hermite,  an  additional  equation  was  introduced  that 

might  require  the  third  derivative.  However,  since  the  first  derivative 

ia  saved  for  the  partition  point,  no  calculation  is  needed  there.  As  a 

result,  the  derivatives  needed  are 

Y"  -  2(,f(x+h)  -  2f(.x)  +  f(x-h))  _  f'(x-Hh)  -  f'(x-h) 
h  "         h2  2h 

vw  =  15(f(x+h)  -  f(x-h))  _  3(f'(x+h)  +  8f'(x)  +  f'(x-h)) 
h-  "       2h3  2h2 

Similarly,   the  trimite   forms   require   the  third  and  fourth  derivatives,  but 

do  not  need  the  first  or  second.      These   are 

Y„,      105(f(x+h)    -f(x-h))       f(f,(x+h)  +  f'^-h))  +I8f'(x) 
8h3 

,    3(f"(x+h)    -  f'(x-h)) 
8h 

Y"h  -  T2(f(x-fh)    -  2f(x)   +  f(x-h))  _  39(f'(x+h)    -  f'(x-h)) 

t  k  ~  3 

h  2ti 

3(f"(x+h)    -  2Uf'(x)   +   f'(x-h)) 
2h2 

In  implementing  the  methods  given  above,  consideration  needs  to  be 

given  to  the  method  of  storage.   In  the  methods  involving  the  generation  of 

an  initial  partition,  it  would  be  possible  to  generate  a  natural  order  of 

points  or  in  some  order  that  is  suggested  by  the  method  of  generating  the 
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points.   If  the  initial  partition  is  generated  in  something  other  than  the 
natural  sequence,  it  might  he  better  to  reorder.   This  question  is  answered 
by  looking  at  the  subroutine  used  to  solve  the  system  of  differential 
equations.   In  this  subprogram  there  are  several  places  where  the  user's 
function  is  needed  and  the  values,  of  the  partition  points  will  be  accessed. 
Each  access  will  require  that  an  order  vector  be  used  to  obtain  the 
positional  value.   There  are  also  several  other  arrays  in  this  subprogram. 
If  these  are  stored  in  the  order  that  the  partition  is  stored,  further  time 
is.  consumed  in  determining  that  address.   If  these  additional  arrays  are 
stored  in  an  order  different  from  that  of  the  partition,  then  some  problems 
will,  occur  whenever  the  partition  needs  to  be  changed.   For  these  reasons, 
it  was  decided  to  maintain  the  natural  order,  regardless  of  how  the  partition 
is  generated. 

2.2  Error  Estimates 


If  the  method  is  going  to  be  automatic,  a  part  of  the  procedure 

must  be  devoted  to  calculating  some  sort  of  estimate  of  the  error  and 

deciding  some  course  of  action  based  on  the  results.   There  are  basically 

three  areas  that  need  to  be  discussed.   Two  of  these  involve  being  able  to 

determine  an  estimate  of  errors  introduced  through  spatial  discretization, 

and  also  an  estimate  of  errors  incurred  by  solving  the  system  of  ordinary 

differential  equations.   The  first  type  of  error  will  be  called  d  and  the 

second  d  .   The  third  area  of  discussion  involves  the  need  to  combine  these 

errors  to  produce  a  result  that  gives  some  measure  of  the  overall  error. 

In  equation  (l),  the  right  hand  side,  A  (u.),  can  be  written  as 

n      1 

f(u. ,   A        (u.),  A  (u. ),    x. ,   t)  (Ik) 

l        n,x     i  n,xx     l  l 

where  A        (u. )    and  A  (u. )    are   analogs   of  the  first   and  second 

n ,x     i  n,xx     i 

derivatives,   respectively.      Putting  this   expression   in  place  of  A   (x. ),    (l) 
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du. 

dT  =  f(V  An,x(V>  An.aCut),  x.  ,  t)  (15) 

The  analogs  A^tu,)  and  A^^fo.)  produce  a  truncation  error 

dx(Xl,t).  An  estimate  of  d^x^tl  can  be  obtained  by  perturbing  the  right 

hand  side  of  equation  (.15)  by  the  error  estimates  of  A   (u  )  and 

n,x  i 

An,xx(ui)  (ei  and   e2'  res-Pec-tively).   This  would  give  the  expression 

f(vAn,x(V  +  er  V(ui}  +  vxi>  *> 

Subtracting  the  expression  (lU)  from  this,  the  result 

dx(xist)  =  f(u.,  Anjx(u.)s  Anjxx(u.)5  x.,  t)  -  f(u.,  An)X(u.)  +er  A^Cu.) 

+  e2>  xi'  ^  (16) 

is  obtained. 

In  order  to  be  able  to  obtain  this  estimate,  the  values  of  e  and 
e2  will  have  to  be  obtained.   If  the  analogs  involved  were  just  those  of  the 
Lagrangian  polynomials,  then  several  approaches  are  available  for  obtaining 
the  error  terms  e±   and  eg.   Since  the  present  approach  is  to  use  Hermite 
type  polynomials  as  well,  the  approach  of  Wendroff  [2]  will  be  followed. 

If  y(x)  is  any  polynomial  approximation  to  f  (x),  then 

f(x)  =  y(x)  +  p(x)  G(x)  (17) 

where  P(x).  G(x)  is  the  error  term,  G(x)  is  an  unknown  function, 

m.       m  m 

P(x)  =  (x-a^  X(x-a2)  2...(x-ak)^: 

and  m  is  the  number  of  values  obtained  at  a  . 
x  i 

In  order  to  find  the  error  estimate  of  f^'(x)  (the  trimite 
method  would  require  that  j=U),  take  the  j-th  derivative  of  both  sides  of  (17) 
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For  the  first  four  derivatives  (assuming  that  G  is  continuous 
and  has  j  continuous  derivatives)  the  results  obtained  would  he 
f'(x)  =  y'(x)  +  P'(x)  G(x)  +  P(x)  G'(x) 

f"(x)  =  y"(x)  +  P"(x)  G(x)  +  2P'(x)  G'(x)  +  P(x)  G"(x)  (iQ) 

f*(x)  =  y"'(x)  +  Pm(x)  G(x)  +  3P"(x)  G'(x)  +  3P'(x)  G"(x)  +  P(x)  G"'(x) 
fH(x)  =  ymt(x)  +  j*(x)  G(x)  +  UF»(x)  G«(x)  +  6P"(x)  G»(x)  +  UP'(x)  G"'(x) 

+  P(x)  G""(x) 

The  first  problem  is  to  determine  the  validity  of  the  assumption 
that  G  G'  G",  G,M,  Gm'  are  actually  continuous  functions.  To  do  this,  for 
G,  solve  (IT)  for  G  giving 

G(x)  =  f(x)  -  yU)  (19) 

Glxi      P(x) 

Since  f(x)  is  assumed  continuous,  G(x)  is  continuous  whenever  P(x)  t   0. 

Consider  now  a  point  a.  ,  p(<3)(a.)  =  0  for  0  <_  J  <  m.  but  P  x    (a.)  t   0. 

It  is  also  noted  that  [f(x)  -  y(x)]  =  0  for  0  <_  J  <  m..   Therefore,  G^) 

is  defined  by 

m.        m. 
f  X(a,)  -  y  1(a.) 


G(  a.  )  =       m 

P  X(a.) 

which  by  L'Hospital's  rule  shows  the  continuity  of  G. 

To  show  that  G' (x)  is  continuous,  Wendroff  [2]  takes  the  derivative 

of  both  sides  of  (19)  giving 

,,  v   P(xUf'(x)  -t'(x)1  -P'(xUf(x)  -y(x)]         (20) 

[pun2 

He  then  defines 

m, 


h(x)-  U^flxf^  (x-a±)  E^ 
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where  k  is  the  number  of  interpolating  points.  Substituting  this  into  (20) 
and  using  L' Hospital's  rule  again,  he  arrives  at  the  continuity  of  G'(x) 

(mi+l) 

provided  f      (x)  is  continuous  for  the  maximum  m  . 

i 

The  same  approach  can  be  used  to  show  that  each  of  G" ,  Gm,  G""  is 
also  continuous  provided  f  has  a  sufficient  number  of  continuous  derivatives. 
Since  the  only  major  difference  is  that  in  each  succeeding  case  the 
computations  become  much  more  involved,  only  the  proof  of  the  continuity  of 
G"  will  be  shown. 

To  this  end,  take  the  derivative  of  both  sides  of  (20)  giving 
_  2(gf)2(f(x)-y(x))-|^f  (f(x)-y(x)) 

u  W  -^       

2(PTxT')  <f,(x)-y'(*))  +  (F"(x)-y"(x)) 

P(x)  (21) 

Again,  the  problem  of  continuity  exists  only  at  points  where 

P(x)  =  0,  since  otherwise  all  functions  are  continuous.   So  consider  a  point 

ai  where  P(x)  has  mi  zeros.   Multiply  the  numerators  and  denominators  of  (21) 

o 
by  (x-a^   giving  (let  f  e  f(x),  y  E  y(x),  P  E  p(x)) 

_  ^*~\)2    (f-)2  (f-y)-(x-a,)2  p(f-y)-(x-a.)2.2.£l  (f'-y>) 

(x-a.)2  P 

(x-a  )2  (f"-y") 

+  =- 


(x-a.)2  P  (22) 

Since  the  denominator  has  m. +2  zeros  at  a.  and  the  (m. +2)-nd 
derivative  is  nonzero,  the  limit  of  G'^a..)  exists— provided  the  numerator 
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vanishes  for  the  first  m.+l  derivatives.   The  continuity  of  G'(ai)  and 

definition  of  the  derivative  then  would  allow  the  conclusion  that  G"(ai) 

is  continuous.   It  is  assumed  that  f(x)  has  a  sufficient  number  of  continuous 

derivatives . 

Define 

pi  k    m 

h(x)  =  (x-a.)-=  (x-a.)   Z  ^j 

and 

d(x)  =  (x-a  )2f-=  (x-a  )2  [   Z    (     ?    )   +     Z        J  \2    1   (23) 
1   P        X  j,£=l  U  ajn   V    J*l  (x-a.T 

It  is  noted  that  h'(x)  and  d'(x)  are  then  given  by 

k     m.     (x-a.  )m. 
j=l  (x-ao}   (x-a/ 


and 


2 
k    2 (x-a.)  m.m,    2 (x-a. )  m.m 

d'(x)  =   Z    [-, 1,      i  \    - i$ ^-H 

j^=l   ^V^V    (x-a,)2(x-a) 

k    m.m.    (x-a.  )  m.m. 
j=l  ^^V      (x-a.)2 

o 
k   2(x-a.  )  m.  (m.-l)    (x-a.)   m.  (m  -1) 
+  z      [ 1    ,1    ,1 _  1 J— J ]      (2U) 

J-l       (X"aJ}  (x-a.)2 

If  (23)  and  (2k)    are  evaluated  at  a. ,  they  become 
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h(a.)   =  nu,   d(ai)   =  m.  (m. -l), 

km  km 

h'(a   )   =     Z     t — J-^r,   d'(aj   =  2m.      Z     7 — J-_ 


Substituting  these  into   (22),  the  numerator  G"(x)  becomes 
2h   (x)    (f-y)   -  d(x)    (f-y)    -  2{x-e,±)  h(x)    (f'-y')   +   (x-a.  }2    (f"-y")        (25) 

The  derivatives   of  each  of  the  four  terms    (call  them  A,   B,    C,    and  D)   in 
(25)  are  now  evaluated  using  Leibnitz'    rule.      This    gives 


1.      A  =   2^-(h2(x)(f-y))    =      Z      0(h2(x))(r)(f-y)(a-r) 
dx  r=0 


This  term  is  now  evaluated  at   a    .      Since    (f-y)   has   m.    zeros,   A  is   identically 


zero   for  0   <   a   <  m 
—  1 


0  (m.  )  9  (m.  ) 

for  a  =  m.  ,   A  =  2«ti(x)(f-y)      X     =  2m? (f-y)      x 

for  a  =  m 


i+l» 


2  (m.+l)  (m.  ) 

A  =  2-h    (x)(f-y)  +   Uh(x)-h'(x)(m.+l  )(f-y)      X 


2           (m.+l)                                       fo)    k  m. 

=  2m? (f-y)    x        +   Urn.  (m.+l  )-(f-y)    X      Z  — J- 

1                                11                             .   ,  x-a . 

j=l  J 


2.      B   5   -^_  (d(x)(f-y))  =     I       (a)(d(x)(r)(f-y)(a-r) 
dx  r=0      r 


B  is  now  evaluated  at   a  .      For  0<a<m. ,   B  =  0 

i  —  1 

for  a  =  m.  , 
1 


(m. )  (m. ) 

B  =  -d(x)(f-y)      x     =   -m.  (m.+l)(f-y)      1 


1      1 


2k 


for  a  =  m.+l, 


(m.+l)  (m ) 

B  =  -(d(x)(f-y)     X         +  (m.+l)   d'(x)(f-y)  ) 

(m.+l)  (m.)    k       m. 

=  -m.  (m.+l)  (f-y)      X  -  (m.+l  )•  2-  m  •  (f-y)   X     Z     t-^—t 

3.      c  =  „2^L  ((x_a.)   h(x)    (f-y'))   =  -2     I      (!)E(x-a.)  h(x)](r)    (f-y)(a+1"r) 
dxa  X  r=0 

Again,    for  0<a<m. ,C=0  since  whenever  r  =  0,   x-ai   =  0. 


for   a   =  m. 

l 


(m  )  (m   ) 

C  =  -2m..h(x)-  (f-y)      X     =  -2m.  (f-y) 

for  a  =  m.+l 

(m+1) 
C  =  -2(m1+l)-h(x)-  (f-y) 

(m.+l)(m   )  (m  ) 

..2._i-^ i-2-h»(x)-(f-y) 

0  (m.+l)  (m  )      k       m 

=  -2(m?+m.)(f-y)  -  2(mf+m,  )    (f-y)  I     -f- 

11  X     X  j=l  X"aj 


and 


d°L  r/„  „    \2/^„  „,ni   -     ?      (*\f(v_0    \2\r < *,Aa+2-T 


k.     D  -S-j  [(x-a.)^(f"-y")]   =     £      ( ")  ( (x-a.)^)r(f-y )' 
Since  D  =  0  when  r  <   2 ,   D  =  0   for  0   <_  a   <  m 


for  a  =  m. 

i 


(m   ) 
D  =  m.(m.+l  )(f-y) 


i     i 
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for  a  =  m.+l 


1 


(m.+l) 
D  =  (m.+l  )-m..  (f-y)   1 

Combining  A,  B,  C,  and  D  for  a  =  m. ,  the  results  are 

r22  ?  P  ^m'  ) 

S  =    [2m.    -  m.    +  m .    -  2mJ   +  mT   -  m   ]    (f-y)      x     =  0 

for  a   =  m.+l 

l 

2  2  2  P  (m.+l  ) 

S  =    [2m     -  m     +  m.    -  2m.    -  2m.    +  mf  +  m   ]    (f-y)      x 


k       m 


2  ?  p  Km 

+    [Urn     +  Urn.    -  2m?   -  2m.    -  2mf   -  2m.  ]    (    E     — <M 
1  x  x  !  i  i        ,  .    x-a/ 

J* 

=   0 

Therefore,  G"(x)  is  continuous. 

Once  the  functions  Gf  ,  G" ,  G«\  and  G""have  been  shown  to  be 
continuous,  it  is  possible  to  obtain  a  value  for  each  of  them.   Since 
the  computations  are  again  very  similar,  it  will  be  shown  for  one  of 
the  four. 

Consider  the  case  of  G»"(x)  .   Define  a  function  H(z)  by 
H(z)  =  (G(x)  -  G(z)  +  G'(x)(z-x)  +  ^G"(x)(z-x)2 

+  \  G«'(x)(z-x)3  +  ^  G»«(x)(z-x)1|)-P(z)  (26) 

Wendroff  [2]  proves  a  generalized  Rolle's  theorem  which  states 
that  if  f(x)  has  an  n-th  derivative  in  [a,b],  n  >_1,  and  if  f(x)  has  a 

zero  of  at  least  m  at  a  ,  where  a  <   a,  :  a0  <...<  a,  <  b,  and  if 

x     1         —  l—  2—   —  k—  ' 


k 


2  m.  ^n+1,  then  there  exists  an  n  e[a,b]  such  that  rn^(n)  =  0. 
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In  the  case  at  hand,  H(z)  has  at  least  n+5  zeros  including  five 
at  z  =  x,  and  n  associated  with  P(z).   This  would  mean  that  there  is  an 
n  such  that  H^    (n)  =  0.  But  on  taking  the  derivatives,  after  multiplying 
out  the  right  side  of  (26)  and  noting  that  P(z)  G(z)  =  f(z)-y(z),  the 
function  "becomes 

H(n+M(n)  =f(*+D(n)  +in±iilGm<x)  =0 

Solving  this  expression  for  G""(x)  gives 

The  general  form  would  then  "be 

o(J)(x)=I^TTf(n+J)(n) 

In  order  to  know  the  value  of  the  error  estimates  given  in  (l8), 
the  values  P,  P'  ,  P",  P"',  P""need  also  to  be  obtained.  Since  the  polynomials 
in  the  present  case  are  the  Lagrangian,  Hermite  and  trimite  polynomials,  P(x) 
takes  on  the  form  p(x),  p2(x),  p  (x)  where  p(x)  =  (x-a^  (x-a2) . .  .(x-an) . 
The  derivatives  of  these  polynomials  are  given  in  Table  1. 
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Lagrangian  Hermite 

P'W  P*(x)  2p(x)p»(x) 

P"(x)  Pn(x)  2[(p'(x))2+p(x)p"(x)] 

F"(x)         not  needed  2[3p'  (x)p"(x)+p(x)p"'(x)  ] 

P"'(x)  not  needed  not  needed 

Trimite 
p'(x)  3p2(x)p'(x) 

P"U)  6p(x)(P'(x))2+3p2(x)p"(x) 

F"(x)  6(p'(x)34-i8P(x)p'(x)p"(x)+3p2(x)p"(x)) 

P"^)  36(p'(x)2p"(x)+l8p(x)(p»,(x))2+2^p(x). 

P'(x)p,"(x)+3(p(x)2p"'(x) 
Table  1.   Expressions  for  Error  Terms 
Considering  that  p(x)  =  0  and  p"(x)  =  0,  since  the  calculations  are  done 
at  the  center  point  of  an  odd  number  of  interpolating  points  (in  this  case 
n=3),  the  substitution  of  the  values  of  Table  1  into  (l8)  yields 
1.   for  Lagrangian  polynomials 

first  derivative  error  estimate  is 

P3U)  f<3),  ,  ,   , 

"51 — f   (n)  (27a) 

second  derivative  error  estimate  is 

2P3U)  Jk).    ,  ,   , 

—IT! — f   (n)  <2Tb) 
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2.   for  Hermite  polynomials 

second  derivative  error  estimate  is 


2(.pi(x))2  ,„ 

— ^ f(6,(n)  (28a) 


third  derivative  error  estimate  is 


6CpUx))2  m 

~ f   (n)  (28b) 


3.   for  trimite  polynomials 

third  derivative  error  estimate  is 

«'^»3  „<9> 


f^'(n)  (29a) 


fourth  derivative  error  estimate  is 

f(l0)(n)  (29b) 


2U^3(x))3  .(10) 


10! 

The  values  of  these  estimates  are  then  put  into  equations  like  (l6)  in 
place  of  e  and  e  to  determine  d  (x. ,t).   In  the  case  of  the  Hermite  or 
trimite,  d  (x  ,t)  would  be  found  by  using  equations  (12)  and  (13)  for 

XI 

the  development  of  d  (x.  ,t). 

In  the  time  direction,  d. (x. ,t)  can  be  found  from  the  automatic 
program  used  to  solve  the  system  of  equations  generated  by  the  discretization, 

Once  these  error  estimates  have  been  determined,  there  are  two 
ways  that  the  overall  error  could  be  approached.   First,  the  sum  of  the  two 
could  be  used  by  asking  that  during  any  time  step  this  sum  would  be  less  than 
the  user's  requested  error.  This  would  be  saying  to  the  user  that  the 
equation  being  solved  is  close  to  his.   However,  the  actual  error  in  the 
solution  will  depend  on  the  nature  of  the  equation.   If,  therefore,  the 
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user  desired  a  particular  error  in  the  solution,  he  would  need  to  have 
some  idea  about  the  nature  of  the  solution  so  that  a  value  of  the  error 
request  can  be  obtained.   This  procedure  would  take  less  computer  time 
but  would  involve  more  thought  for  the  user. 

The  second  approach  would  be  to  attempt  to  obtain  an  estimate 
of  the  error  in  the  solution.   This  could  be  done  by  numerically  solving 
a  system  of  asymptotic  error  equations  along  with  the  discrete  system  of 
differential  equations.   The  system  of  equations  being  solved  is  of 
the  form 

du. 


The  system  differs  from  the  desired  system  in  that  the  errors  d  and 

dx  have  been  introduced.   This  would  lead  to  a  system 

3u(x.  ,t) 

~ =f(u(x.,t),  AnfX(u(x.,t)),AnjXx(u(x.,t)),  x.,  t)  +  dt  +  dx 

Subtracting  and  letting  e.(t)  =  u.(t)  -  u(x. ,t),  the  equations  obtained  are 

de. 

dT^V  An,x(ui}'  An,xx(ui}'  V  V   ei  "  dt  - d   (30) 

Since  the  solution  to  this  equation  need  not  be  accurate,  it  may  be  possible 

to  use  a  lower  order  method  to  solve  this  last  system  provided  f   (u(t))  is 

li- 
the same  for  both  systems. 

The  advantages  of  the  latter  method  led  to  solving  the  system 
(30)  along  with  the  original  system. 
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2.3  Efficiency 

Once  methods  are  available  for  solving  the  problem  and  for 
estimating  the  errors,  attention  needs  to  be  given  to  producing  an 
efficient  system.   That  is,  the  solution  needs  to  be  found  in  such  a  way 
that  it  has  certain  error  properties  and  such  that  the  work  involved  in 
obtaining  that  solution  satisfies  some  sort  of  minimization  of  computer 
time.   If  greater  accuracy  is  desired,  it  is  expected  that  more  computer 
time  will  be  used.   If  there  is  a  decrease  in  computer  time,  then  provided 
the  overhead  is  always  kept  at  a  minimum,  a  loss  of  accuracy  is  expected. 

The  goal  would  be  to  have  the  error  in  the  result  less  than, 
but  as  close  as  possible  to,  the  amount  requested  by  the  user  in  such  a 
way  that  no  more  computer  time  is  used  than  is  absolutely  necessary. 
The  primary  emphasis  needs  to  be  that  of  keeping  any  estimated  error  less 
than  the  user's  request. 

If  for  some  reason  the  estimated  error  becomes  too  large ,  over- 
head will  have  to  be  expended  to  rectify  the  problem.   If  the  error 
becomes  much  smaller  than  the  user's  request,  it  might  be  desirable  to 
modify  the  method — provided  the  overhead  is  less  than  the  expected  savings 
in  time . 

There  are  three  ways  to  satisfy  these  goals: 

1.  modify  d  (x. ,t)  by  changing  the  error  request  given  to  the 
subprogram  solving  the  system  of  equations  (DIFSUB), 

2.  modify  d  (x. ,t)  by  changing  the  spacing  in  the  partition, 

3.  modify  d  (x. ,t)  by  changing  the  order  of  the  spatial  analog. 

.A.    J- 
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When  "beginning  the  integration,  it  would  be  desirable  to 
determine  first  (if  possible)  those  items  involving  the  greatest  amount 
of  overhead.   After  the  integration  has  proceeded  a  number  of  steps,  the 
item  involving  the  least  amount  of  overhead  or  the  item  least  apt  to 
adversely  affect  the  solution  should  be  changed  first.   Care  must  also  be 
taken  that  changes  do  not  occur  too  frequently  since  modifications  will, 
in  general,  have  a  tendency  to  introduce  additional  errors. 

Of  the  three  items  above,  the  first  has  the  most  desirable 
properties.   The  change  of  this  parameter  could  cause  the  vectors  in 
DIFSUB  to  be  multiplied  by  some  factor  depending  on  the  change  required 
in  step  size.   On  the  other  hand,  the  last  two  involve  a  greater  amount 
of  overhead  and  also  involve  interpolation  that  could  introduce  larger 
errors  than  those  involved  with  changing  the  time  parameter. 

Consider  the  problem  of  obtaining  an  initial  partition  or  of 
changing  to  a  new  partition.   The  problem  is  to  determine  the  order  of  the 
method  and  the  points  to  include  in  the  partition  so  that  the  integration 
may  be  accomplished  most  efficiently  and  yet  have  the  desired  accuracy. 
If  it  is  assumed  that  given  a  partition,  the  majority  of  work  is  involved 
with  the  function  evaluation  or  calculation  of  derivatives ,  and  in  the 
solution  of  a  system  of  equations  required  by  DIFSUB,  it  would  be 
possible  to  develop  a  sort  of  work  estimate  based  on  the  number  of 
multiplications  and  divisions  present. 

In  the  case  of  the  derivative  calculation,  the  two  things  to 
consider  are  the  calculation  of  the  derivative  and  the  calculation  of  its 
error  estimate.   In  the  first  case,  we  are  using  the  Lagrangian,  Hermite, 
and  trimite  three-point  formulas,  so  that  the  coefficients  are  known  as  socn  as 
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the  step  size  is  given.   In  the  case  of  the  error  estimate,  the  coefficients 
for  its  interpolating  polynomial  are  calculated  and  then  stored  for  a  given 
partition.  When  this  estimate  is.  calculated,  there  -will  he  four  inter- 
polating points  for  each  derivative  (-with  the  exception  of  the  second 
derivative  using  the  Lagrangian  method  which  requires  five  points).   The 
following  table  lists  the  number  of  operations  used  in  calculating  the 
derivatives  for  each  method.   The  derivatives  requiring  calculation  are 
the  first  and  second  for  the  Lagrangian,  second  and  third  for  Hermite,  and 
third  and  fourth  for  the  trimite. 


Multiplications/Divisions 
In  Derivative 


Multiplications/Divisions 
In  Error  Estimate 


Lagrangian 

Hermite 

trimite 


lower  higher 

h  5 

8  8 

12  12 


lower  higher 
1      2 
k  5 

7     8 
Table  2 .   Number  of  Operations  in  Equation  Evaluation 

To  aid  in  determining  which  method  would  be  most  efficient,  the 
user  is  to  indicate  the  derivatives  required  for  each  of  the  methods  used. 
In  solving  the  system  of  equations,  the  procedure  is  divided 

into  two  steps : 

1.  decomposing  the  matrix  into  upper  and  lower 

triangular  matrices, 

2.  using  the  decomposed  matrix  to  solve  the  system  of  equations. 
In  each  of  these  steps ,  the  number  of  nonzero  multiplicative  operations 
required  for  each  equation  is  proportional  to  one-half  the  number  of 

of f -diagonals .   In  the  case  of  decomposing  the  matrix,  this  proportionality 


33 

constant  is  two.   In  the  case  of  the  solving  step,  the  constant  is  three. 
The  number  of  off -diagonals  is  two,  five,  and  eight  for  the  Lagrangian, 
Hermite,  and  trimite,  respectively.   The  average  number  of  operations  per 
equation  for  each  method  with  each  step  is  given  in  Table  3. 

Decomposition  Solving 

Lagrangian  2  o 

Hermite  5 

trimite  8  12 

Table  3.   Number  of  Operations  in  Decomposing  and  Solving 

The  problem  in  trying  to  bring  the  three  sets  of  numbers  together 
is  that  the  relative  frequency  of  the  number  of  function  (Nf )  ,  decomposition 
(Nd),  and  solving  (Hq)  calls  depends  on  the  nature  of  the  problem.   The 

number  of  function  evaluations  can  be  given  directly  as  a  function  of  N 

d 

and  N  bv 

s  J 

N  =  n  *N  +  i.  N 
f    m  d   2   s 

where  r^  is  the  number  of  equations  in  the  system.   The  reason  for  the 

factor  of  \  is  that  for  each  call  to  obtain  the  function  evaluation,  the 

solving  routine  is  called  for  use  by  the  system  of  differential  equations 

and  for  use  by  the  system  of  asymptotic  error  equations.   The  first  term 

arises  from  the  fact  that  DIFSUB  calculates  a  Jacobian  by  a  numerical 

differencing  method  that  calls  for  the  function  evaluation  for  each  equation 

in  the  system.   The  ratio  (R)  of  Nd  to  Ng  is  variable  and  depends  on  the 

solution  to  the  system  of  equations.   In  systems  where  the  Jacobian  varies 

slowly,  the  ratio  ranged  from  about  .033  to  .0^5.   Since  those  equations 

were  very  well  behaved  and  since  an  increased  ratio  indicates  an  increased 
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amount  of  work,  a  value  of  .05  will  he  used  initially.   The  ratio  will 

then  he  modified  as  the  integration  progresses .   At  the  start  of  the 

integration,  the  value  N,  could  he  written  as 

N,  =  .05  N 
d        s 

An  estimate  of  the  work  required  for  each  method  can  now  he 

obtained  as  a  function  of  the  number  of  calls  to  the  solving  routine. 

Such  an  estimate  is  given  by 

W  =  n  {^7-W..  +  R(W  •-—  +  W.)  +  W  }«N 
m    m  2M  fl      f2  M     d     s   s 

where  n  is  the  number  of  equations,  W,  and  W  are  the  number  of  operations 
m  as. 

in  the  decomposition  and  solving  subroutines.   M  is  the  number  of  values  stoi 
at  each  point,  and  W   ,  W   are  the  operations  required  in  a  function  call  fc 
a  step  in  the  solution,  and  a  function  call  in  determining  the  Jacohian.   The 
latter  distinction  is  brought  about  because  the  process  of  calculating  the 
Jacobian  does  not  require  any  error  estimate  calculations. 

In  order  to  decide  which  method  should  he  used,  the  values  of 

W  /N  for  each  of  the  methods  are  compared.   Although  one  is  interested  in 
m  s 

using  the  method  with  the  smallest  value  of  W  /N  ,  there  is  an  additional 

m  s 

complication  in  trying  to  decide  which  method  to  use.   This  involves  the 
fact  that  if  an  attempt  were  made  to  change  the  partition  during  the 
integration  to  a  higher  ordered  method,  the  interpolation  polynomials  would 
be  such  that  error  estimates  would  be  too  large  to  work  with.   As  a  result, 
it  might  be  appropriate  to  bias  the  selection  slightly  in  favor  of  the  highei 
ordered  methods  during  the  initial  partitions  selection.   After  integration 
has  begun,  if  a  change  is  to  be  made  in  the  partition,  the  simpler  lower 
ordered  methods  would  be  preferred. 
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Once  the  partition  has  been  selected,  it  would  then  be  possible 
to  adjust  dt  to  some  extent  so  as  to  keep  the  error  close  to  the  user's 
request.   If  DIFSUB  is  asked  to  produce  a  large  dt  then  it  will  increase 
the  step  size  or  lower  the  order  of  integration.   If  asked  to  produce  a 
smaller  error,  it  will  decrease  the  step  size. 

The  modification  of  dt  and  dx  should  not  occur  too  frequently 
since  the  introduced  errors  could  become  large  enough  that  a  satisfactory 
solution  is  impossible. 


2.k     Initial  Partiti 


on 


In  order  to  start  the  integration,  it  is  necessary  to  have  an 
initial  partition.   Since  any  modification  in  the  partition  will  introduce 
additional  errors,  it  is  important  to  find  a  partition  that  will  be  the 
basis  of  an  efficiently  obtained  solution  that  will  satisfy  the  user's 
request.   One  way  to  develop  such  a  partition  would  be  to  control  the 
amount  of  error  present  in  the  derivative.   This  would  lead  to  trying  to 
obtain  a  spatial  error  parameter  that  would  be  used  to  calculate  this 
initial  partition. 

There  are  three  basic  approaches  that  could  be  followed.   One, 
the  user  could  be  asked  to  provide  this  parameter.   Two,  a  set  function  of 
the  user's  overall  error  request  could  be  used  to  develop  such  a  parameter. 
Three,  an  arbitrary  set  of  partitions  using  equal  spacings  could  be  set 
up  and  then  tested  to  see  what  error  results  in  the  asymptotic  error  formula 

The  first  method  of  the  three  has  the  least  desirable  properties 
from  the  standpoint  of  the  user.   It  requires  that  he  make  some  sort  of 
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judgment  about  the  equation.   If  his  judgment  gives  poor  results,  then  the 
parameter  would  have  to  be  modified.   Since  the  program  would  eventually 
have  to  modify  this  parameter,  it  would  be  appropriate  to  generate  the 
parameter  initially. 

The  second  method  has.  the  advantage  that  the  function  could  be 
an  easily  calculated  value.  For  example,  one  could  have  this  parameter 
be  one-tenth  the  value  of  the  user's  request.   If  this  generates  a 
partition  that  creates,  valid  errors,  then  no  unnecessary  work  has  been 
done .   If  the  errors  generated  are  too  large ,  then  the  partition  would 
have  to  be  recalculated  based  on  some  formula  of  parameter  modification. 

The  third  method  would  start  with  a  partition  using  a  step  size 

a— b 

of  ,  where  a  and  b  are  the  end  points.   Using  this  partition,  the 

2m 

integration  would  be  performed  over  one  step.   The  error  in  the  solution 
would  then  be  compared  with  the  errors  in  the  derivatives.   The  process 
would  be  continued  until  the  errors  in  the  one  step  of  integration  correspond 
with  the  user's  request.   The  maximum  error  estimate  in  the  derivatives 
would  then  be  used  to  generate  the  working  partition.   The  disadvantage  with 
this  approach  is  that  overhead  will  always  exist.   If  the  first  partition 
yields  a  value  for  the  spatial  error  parameter,  the  working  partition  would 
still  need  to  be  generated.   The  exception  to  this  would  be  the  case  where 
the  working  partition  contains  equal  spacings. 

Since  the  second  method  is  simpler  to  work  with  and  since  it  can, 
if  required,  perform  somewhat  like  the  third  method,  it  was  decided  that 
the  second  approach  would  be  used. 
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Once  an  error  parameter  has  been  determined,  the  working  partition 
needs  to  be  generated.   Two  approaches  were  considered.   One,  it  would  be 
possible  to  start  at  one  of  the  end  points  and  then  choose  each  successive 
point  so  that  the  error  condition  is  met,   Two,  it  would  be  possible  to 
start  with  three  intervals  and  then  subdivide  each  interval  in  which  the 
error  condition  is  not  met. 

In  the  first  method,  one  would  start  by  adding  two  points  to  the 
partition,  a  center  point  and  a  second  end  point  of  the  first  subinterval. 
Thes.e  two  points,  would  be  chosen  so  that  the  error  condition  would  be  met. 
The  next  step  would  be  to  use  the  second  point  above  as  the  center  of  a  new 
subinterval,  and  if  necessary,  add  points  on  either  side  of  this  point  so 
that  the  symmetric  three-point  formula  would  again  satisfy  the  error 
condition.   The  process  is  continued  until  no  new  subintervals  need  to  be 
added  to  completely  cover  the  entire  interval  [a,b]. 

Some  of  the  disadvantages  of  this  procedure  are 

1.  each  time  a  new  center  point  is  added,  it  may  be  necessary 
to  backtrack  in  the  interval  to  handle  newly  added  points, 

2.  even  if  all  step  sizes  used  in  the  subinterval  are  powers  of 
two,  there  may  develop  problems  terminating  at  the  second 
endpoint , 

3.  the  possible  backtracking  could  create  programming  complications 

In  the  second  method,  start  with  the  three  points  P  =  a  +  ^^  P  =  n  +  ^ZiL 

1        2  '   2   a         k    ' 
b-a 
p3  =  a  +  3'  -jp.   (see  Figure  l) 


a 

P 

I 

I 
2 

P5 
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1                        F 
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Figure  1 
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The  subinterval  [a,P  ]  is  next  checked  to  see  if  the  error  is  reasonable. 

b-a 
If  not,  points  are  added  on  either  side  of  P2  with  step  size  -g— .  If  the 

error  was  acceptable,  then  this  fact  is  recorded  so  that  on  future  passes 

through  the  interval,  this  subinterval  will  be  bypassed.  The  next  step 

would  be  to  check  the  subinterval  [P0,P  ]  to  see  if  the  analog  for  P1 

satisfies  the  error  condition.  At  the  same  time,  the  newly  added  points 

P.  and  Pq  are  checked  to  make  sure  that  they  also  satisfy  the  error 

condition.  This  process  is  repeated  until  all  subintervals  satisfy  the 

error  condition . 

In  order  to  find  out  if  a  given  derivative  satisfies  the  error 
requirements,  it  is  necessary  to  generate  some  sort  of  error  estimate.  An 
attempt  is  made  to  duplicate  the  type  of  error  estimate  that  would  be  made 
in  the  actual  integration. 

In  an  attempt  to  guard  against  unusual  initial  conditions,  the 

b— a   .,       . . , 
testing  at  each  point  is  done  twice — once  with  an  h  =  — —  and  once  with 

2 

h  =  a— —-where  .9  <  a  <  1. 
2n 

2.5.   Changes  in  Integration 

Once  the  integration  parameters  have  been  set,  nothing  needs  to  be 
changed  until  one  of  two  conditions  occur.   If  the  estimated  error  ever  becom 
larger  than  the  desired  value,  then  some  sort  of  corrective  action  will  be 
required  immediately.   If  the  estimated  error  becomes  much  smaller  than  the 
requested  value,  it  may  indicate  that  a  change  is  desirable. 

Consider  the  case  where  the  estimated  error  is  larger  than  desired. 
There  are  two  possible  actions — either  the  error  parameter  to  DIFSUB  or 
partition  could  be  changed. 
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The  first  action  to  be  taken  would  "be  to  attempt  to  change  d  . 
This  would  be  done  provided  the  difference  between  the  estimated  error  and 
the  error  request  is  smaller  than  the  error  reduction  possible  from  DIFSUB 
by  changing  d^.   If  this  condition  cannot  be  met,  or  if  several  changes  in 
this  parameter  have  already  occurred,  then  it  will  be  necessary  to  change 
the  partition. 

Although  it  would  be  possible  to  add  selected  points,  it  was 
felt  that  since  the  changing  of  the  partition  is  being  done  as  a  last 
resort,  the  desirable  action  would  be  to  change  the  entire  partition  to 
one  with  smaller  step  sizes.   There  are  two  other  advantages  in  this 
course  of  action.   First,  the  program  used  to  change  the  partition  could 
be  basically  the  same  as  that  needed  to  modify  it.   Second,  since  the 
entire  partition  now  is  satisfying  a  smaller  value  for  d  ,  there  is  the 
possibility  that  the  partition  would  not  need  to  be  changed  as  frequently. 

If  the  conditions  above  occur  during  the  first  several  steps,  then 
it  is  assumed  that  a  finer  partition  should  have  been  used  initially.   Rather 
than  create  a  new  partition,  it  was  felt  that  a  better  solution  would  be 
obtained  if  the  process  were  restarted  with  a  smaller  value  in  the  error 
parameter.   If  this  process  of  restarting  must  be  repeated  a  few  times, 
then  it  is  assumed  that  the  problem  cannot  be  solved  with  this  method. 

If  the  estimated  error  is  much  smaller  than  that  required,  then 
a  similar  sequence  of  steps  would  occur  as  those  just  discussed.   The  first 
step  would  be  to  see  if  a  change  in  the  time  error  parameter  is  appropriate. 
The  maximum  size  allowed  for  the  time  error  parameter  would  be  approximately 
nine-tenths  of  the  user's  request.   If  the  error  estimate  continues  to 
remain  small,  then  consideration  will  be  given  to  changing  the  partition. 


ko 

Since  going  to  new  partitions  can  introduce  errors,  there  is  an  attempt 

to  delay  any  changes  to  coarser  partitions  unless  there  is  ample  indication 

that  such  a  partition  will  not  have  to  be  changed  back. 

Just  as  in  going  to  finer  partitions,  it  would  be  possible  to  set 
up  the  program  to  delete  selected  points  where  the  error  estimates  indicate 
it  could  be  done.  However,  since  precautions  are  taken  in  not  changing  too 
soon,  it  is  felt  that  when  a  change  is  made  to  a  coarser  partition,  the 
entire  partition  should  be  changed.   Also,  it  would  be  appropriate  to  see 
if  a  lower  ordered  method  could  be  used.   The  ideal  situation  in  trying  to 
go  to  a  lower  ordered  method  would  be  to  be  able  to  keep  the  same  partition 
but  just  reduce  the  number  of  items  saved  at  each  point.  Normally,  when 
changing  to  a  lower  ordered  method,  one  would  expect  to  have  to  interpolate 
some  in  order  to  make  up  for  the  loss  of  values  at  the  partition  points. 

2.6  The  User's  Equations 

The  communication  between  the  user  and  this  automatic  system  of 
subprograms  used  to  solve  the  user's  equation  (called  the  system)  should 
be  kept  as  simple  as  possible.  Besides  a  program  to  drive  the  system,  the 
user  would  be  expected  to  provide  subprograms  that  give  the  differential 
equation,  initial  values,  and  boundary  values. 

In  the  case  of  the  differential  equation,  the  user  would  be 

expected  to  provide  an  equation  of  the  form 

2 
3u    „/    3u   8  u     .  v 

W=  f(u>  3?  72'  X'  t}' 
9x 

The  system  would  provide  the  values  for  the  arguments  of  the  function  f ,  and 

his  program  would  return  the  time  derivative. 
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If  the  higher  ordered  methods  are  also  used,  then  two  other 
equations  would  also  be  included.   These  equations  would  have  the  forms 

9v  _   ,      9v  92v     ^ 
9x 

and 

3w  _  ,  ,         9w  92w      , 

9t~~  n(u»  V>  W>  ii"»  — *  x>  t) 

9x 


where 


and 


9x 


-  iJi      92u 

V  "  9x'  V  '"   2 
9x 


Again,  the  system  would  provide  the  values  for  the  variables  of  each 

function,  and  the  user  would  return  the  time  derivatives. 

A  similar  approach  is  used  for  the  initial  values.   Considering 

that  all  ordered  methods  might  be  used,  initial  values  would  be  required 

for  u,  v,  and  w.   If  u(x,0)  ■  r(x),  then  v(x,0)  =  s(x)  =   dr^x)  and 

dx 
2 

w(x,0)  =  t(x)  =  — ^ —  =  —^ — .   If  the  derivatives  of  one  of  these  do  not 
exist,  then  either  an  approximation  with  continuous  derivatives  will  need 
to  be  developed  or  the  method  will  have  to  be  restricted  to  those  lower 
ordered  methods  for  which  the  initial  values  exist. 

The  problem  of  the  boundary  conditions  is  more  difficult  to  handle 
The  general  form  of  the  boundary  condition  is 

u(0,t)  +  u  (0,t)  -  C(t)  =0  at  x  =  0 
and 

u(l,t)  +  ux(l,t)  -  d(t)  =0  at  x  =  1 
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The  user  supplies  these  in  terms  of  algebraic  equations  of  the  form 

r  =  u.  +  A(u. ,  x.  ,  t)  -  d(t) 
1      1   i 

where  A(u.   x. ,  t)  is  some  representation  of  ux  and  r  will  he  the  error 
in  satisfying  the  boundary  condition. 

If  the  Lagrangian  method  is  used,  the  function  A  will  have  to  be 
determined  by  some  form  of  interpolation.   In  the  case  of  the  two  higher 
ordered  methods,  u  =  v.  so  that  A(u±,  x±,   t)  =  v±   could  be  used. 

Since  the  higher  ordered  methods  require  additional  values  at  the 

boundaries,  it  will  be  necessary  for  the  user  to  provide  additional 

equations  involving  these  variables.   It  would  be  ideal  if  these  boundary 

value  equations  could  be  obtained  from  the  boundary  condition  of  the 

original  equation.   The  simple  example  of 

u  =  u  ,  u(x.O)  =  sin(nx),  u(0,t)  =  u(l,t)  =  0  . 
t    xx 

shows  that  this  will  not  work.  The  solution  of  this  equation  is  given  by 

2 

u(x,t)   =  e_7r       sin(   x)  .      From  the   definition  of  v(x,t),  the  solution  is 

found  to  be 

v(x,t)   =   TTe"  cos(ttx) 


Thus, 


-ir2t 
v(0,t)   =  -v(l,t)   =  Tre 


which  cannot  be  obtained  from 

u(0,t)   =  u(l,t)   =  0 
One  possibility  for  handling  the  problem  would  be  to  use  an 
equation  of  the  form 


r  =  v.    -  B(u. ,   v. ,   x. ,   t) 
l  i'      l        l 


where  B  is  an  approximation  of  v.  which  is  developed  from  interior  points. 
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3.   PROGRAMS 

The  purpose  of  this  section  is  to  provide  a  description  of  the 
programs  used  in  developing  the  automatic  solution.  There  are  eleven  such 
programs.   Eight  of  these  programs  (PARDIF,  XSTP,  DIFFUN,  DERIV,  DIVERR, 
DEC0MP,  S0LVE,  DIFSUB)  are  set  programs,  that  make  up  the  method.   The 
remaining  three  (BDRY,  INIT,  UFUN)  are  user-supplied  and  provide 
information  about  the  user's  equation,  initial  values,  and  boundary 
conditions .   The  user  also  needs  a  driving  program  which  would  handle 
the  printing. 

In  order  to  use  the  system,  the  driving  program  calls  PARDIF 
for  each  step  in  the  time  direction.  As  a  result,  the  user  can  determine 
when  and  what  information  he  desires  to  print  out,  and  when  to  stop 
the  integration.   These  programs  are  in  Appendix  A. 

3.1  PARDIF 

PARDIF  is  the  main  program  of  the  method.   It  has  the  basic 
purpose  of  integrating  one  step  for  each  call  made  to  it.   In  performing 
this  step,  it  must  cause  an  initial  partition  to  be  generated,  create 
functional  values  for  this  partition,  integrate  the  step,  decide  whether 
the  step  was  acceptable,  and  decide  whether  modifications  should  be  made. 

The  creation  of  a  partition  is  done  in  two  steps.   First,  it 
calls  XSTP  which  provides  as  output  a  pointer  to  indicate  the  proper  order 
of  points,  and  a  vector  which  contains  partition  values  as  well  as  an 
indication  of  the  other  points  which  are  used  in  analog.   The  second  step 
consists  of  placing  values  in  the  proper  order  in  the  array  x  and  placing 
two  values  in  a  second  array  (INT).   The  one  value  indicates  the  relative 
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position  of  one  of  the  other  points  in  the  analog.  This  will  provide 
sufficient  information  since  the  partition  is  generated  in  such  a  way 
that  two  of  the  three  points  of  the  analog  are  adjacent.  The  pointer, 
therefore,  is  used  to  indicate  the  third  point.  The  other  value  in  INT 
is  the  starting  position  for  the  interpolating  polynomial  of  the 
error  estimate. 

Once  the  partition  is  well-defined,  the  value  of  the  function 
at  each  of  these  points  needs  to  be  found.   If  this  is  the  initial 
partition,  then  the  values  are  obtained  from  an  initial  value  subprogram 
(INIT).   If  the  partition  is  a  modification  of  an  existing  one,  the  program 
tries  to  find  the  values  in  the  old  partition.   If  this  attempt  is 
unsuccessful,  then  the  values  are  found  by  interpolation  by  calling  a 
subprogram  (DERIV) .  . Since  there  are  several  different  vector  components 
that  would  use  the  same  interpolation,  DERIV  is  asked  to  return  the 
coefficients  of  this  polynomial.  The  interpolation  of  all  components  is 
then  accomplished  in  PARDIF. 

Every  time  a  new  partition  is  generated,  the  program  calculates — 
by  calling  DERIV — the  coefficients  used  to  calculate  the  error  estimates. 
These  coefficients  are  then  stored  in  an  array  (DIVER). 

If  all  the  values  associated  with  the  partition  have  been  defined, 
or  if  this  step  does  not  modify  any  of  the  parameters,  the  program  is  ready 
to  integrate  and  check  the  results.   The  integration  is  done  by  calling 
DIFSUB. 

On  returning  from  DIFSUB,  the  program  first  checks  the  return 
code  to  find  out  if  the  step  was  completed  satisfactorily.   If  not,  the 
program  first  tries  to  decrease  d  by  decreasing  the  time  based  error 
parameter  (TPS).   If  this  procedure  is  unsuccessful,  an  attempt  is  made  to 
change  the  partition. 
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Whenever  the  return  code  from  DIFSUB  indicates  that  the  integration 
was  successful,  PARDIF  determines  whether  any  additional  action  is  to  he 
taken.   To  accomplish  this,  the  program  determines  the  maximum  error  (RMAX) 
introduced  in  the  step  by  looking  at  the  results  from  the  asymptotic 
error  equation. 

First,  consider  the  case  where  RMAX  is  larger  than  the  user's 
request  (EPS).   First  of  all,  the  number  of  steps  is  checked  to  see  if  it 
is  within  ten  steps  of  t  =  0.   If  so,  the  time  based  error  parameter  (TPS), 
and  the  space  based  error  parameter  (SPS)  are  reduced  by  a  factor  of  ten, 
and  the  integration  is  restarted  from  t  =  0.   This  process  is  repeated 
three  times  before  taking  a  terminal  exit. 

When  RMAX  is  larger  than  EPS,  an  attempt  is  made  to  decrease  TPS 
until  it  has  a  value  less  than  10~5.   Once  TPS  has  reached  this  value,  the 
space  parameter  (SPS)  is  decreased  by  a  factor  of  . 5  in  an  attempt  to  change 
the  partition.   The  desire  is  to  add  enough  points  to  cause  the  error 
estimate  to  be  less  than  EPS. 

In  the  event  that  RMAX  is  less  than  EPS,  consideration  can  be 
given  to  modifying  parameters  so  as  to  decrease  the  work.   The  program 
DIFSUB  will  automatically  change  the  order  of  the  method  it  uses  and  the 
step  size  associated  with  that  method  to  keep  the  error  close  to  parameter  TPS. 

No  changes  will  be  attempted  in  TPS  or  SPS  unless  two  conditions 
are  satisfied.   One,  the  last  modification  has  not  occurred  in  the  previous 
ten  steps.   Two,  the  error  introduced  in  the  step  must  be  less  than  some 
ot  times  EPS. 

If  the  two  conditions  have  been  met,  TPS  is  increased  by  the 
minimum  of 
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logTPS  +   lofiEPS  gmd  Tps  +  i.  (EPS.RMAX) . 
2  2 

Although  many  other  methods  of  changing  are  possible,  this  choice  appears 

to  change  TPS  rapidly  without  causing  the  errors  to  become  too  large. 

After  three  such  changes  are  made  in  TPS,  consideration  is  given 

to  changing  the  partition.   If  SPS  is  changed  (and,  therefore,  the  partition), 

an  additional  restriction  is  imposed.   This  restriction  is  that  RMAX  be 

less  than  a2 -EPS.   If  this  condition  is  met,  then  SPS  is  given  the  value  of 

the  minimum  of  2-SPS,  SPS  +  ,5(EPS-RMAX)  and  .U-EPS.   If  a  change  does  take 

place  in  SPS,  then  TPS  will  be  given  a  new  value  by 

TPS  =  . 5 • ( EPS-SPS ) . 

Again,  the  changes  represent  one  choice  of  many  possibilities.  This  time 

it  is  desirable  to  increase  the  parameter  rather  slowly  so  that  the  partition 

does  not  change  too  much. 

Consideration  will  also  be  given  to  changing  the  partition  without 
changing  SPS  provided  the  following  set  of  conditions  is  satisfied, 

1.  at  least  30  steps  have  occurred  since  the  previous  modification, 

2.  RMAX  has  remained  less  than  EPS, 

3.  each  component  of  the  spatial  truncation  error  is  less 
than  a2 -EPS. 

3.2  XSTP 

XSTP  has  the  function  of  generating  the  positional  values  of  the 
partition  under  two  different  situations .   The  first  occurs  when  the  initial 
partition  is  generated  and  the  second  whenever  an  existing  partition  is  modif. 
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In  both  cases,  the  program  begins  with  partition  points  at  .5, 

.25,  .75,  and  a  step  size  h  =  £.   On  each  succeeding  pass  through  the 

interval,  h  is  divided  by  2.   Until  h  <_  HMAX  (the  maximum  step  size  the 

user  desires),  the  program  places  points  at  (2i-l)«h,  1  <  i  <  — .  After 

—   —  2h" 

h  >  HMAX,  the  program  places  points  only  when  the  error  estimate  in  the 
analog  does  not  satisfy  the  error  condition.   These  additional  points  are 
placed  on  either  side  of  the  test  point.  Because  of  the  way  the  partition 
is  generated,  each  point  that  is  added  vd.ll  be  halfway  between  two 
existing  points . 

If  the  analog  at  a  particular  point  satisfies  the  error  request, 
this  fact  is  recorded  in  the  vector  NCK.  Before  trying  to  add  points,  the 
program  first  checks  NCK  to  see  if  the  testing  can  be  ignored.   When  all 
values  of  NCK  are  set  to  +1,  then  the  procedure  is  terminated. 

The  differences  between  the  two  uses  of  XSTP  take  place  in 
obtaining  the  functional  values  and  the  error  estimates.   In  obtaining  the 
functional  values  for  an  initial  partition,  the  program  calls  the  user's 
initial  value  program  (INIT).   The  case  of  modifying  an  existing  partition 
is  accomplished  by  the  interpolation  of  values  of  the  previous  partition. 

For  estimating  the  error  of  an  initial  partition,  the  program  adds 
two  points  whose  positions  relative  to  the  points  of  the  analog,  are 
fixed.   Tliis  allow-,  the  coefficients  to  be  calculated  for  all  cases.   When 
the  partition  is  changed,  the  functional  values  are  not  as  easily  determined, 
so  the  error  estimate  is  obtained  by  using  values  determined  for  the  partition 
and  calling  the  subprogram  DERIV  to  interpolate. 

Once  the  partition  for  a  particular  order  has  been  developed,  the 
work  estimate  is  obtained  using  the  ideas  of  section  2.3.   The  entire 
process  is  then  repeated  for  the  next  ordered  method. 
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When  the  partitions  for  all  possible .orders  have  been  developed, 

the  work  estimates  are  compared  to  determine  which  method  should  be  used. 

This  is  done  by  asking  that  between  two  adjacent  ordered  methods,  the 

smaller  of  aw    w   (where  m  is  the  lower  order)  indicates  the  order 
m].   m2         1 

preferred.  The  a  is  used  to  set  the  biases  desired  in  the  selection 

of  a  method. 

Once  the  method  has  been  chosen,  the  program  then  checks  each 
analog  to  make  sure  that  the  smallest  available  step  size  is  being 
used.   It  also  prepares  output  vectors  indicating  the  natural  order  of 
points  and  the  points  that  will  be  used  in  each  analog. 

3.3  DERIV 

DERIV  is  a  double  purpose  subprogram  involved  with  interpolation. 
The  one  task  is  to  calculate  the  coefficients  for  the  error  estimation 
polynomial  of  a  given  partition.  The  other  task  is  to  provide  interpolation 
polynomials  for  calculating  new  values  of  a  partition. 

The  error  estimates  needed  are  given  by  (27),  (28),  and  (29).  Each 
of  these  involves  the  factor  P'(x),  which  can  be  determined  as  a  function  of 
the  spacing  h.   Since  n  =  3,  P  (x)  =  (x-a^  (x-a2)  (x-a3) ,  and  x  =  &2,   the 
value  of  P' (x)  is  given  by 

P'(x)  =  (x-a  )(x-ag)  +  (x-a1)(x-a3)  +  (x-a2)(x-a3) 
=  (x-a  )( x-a _) 

=  -h2 
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Thus,  the  error  terms  become 


-T-  f        (n)    and  -rr  f        (n.)   for  the  Lagrangian 
|gg-  f^(n)   and  1^-  rT)(n)   for  the  Hermite 


and 


ZiMo   f(9)(ri)  md  151200  f(l°}^)  for  the  trimite 
The  problem  remaining  in  determining  the  estimate  is  the 
evaluation  of  f   (n).   In  "the  case  of  the  Lagrangian  polynomials 

f(j)(x)  =  E  *fj)(x)  f(x.)  (31) 

i=l  X        X 

In  order  to  obtain  an  estimate  of  f   (x) ,  and  f   (x),  the  &.(x)'s 

were  chosen  so  that  the  proper  derivatives  existed.   In  the  case  of  the 

third  and  fourth  derivatives,  k   and  5  points  were  required,  respectively. 

Using  these  polynomials  gives 

f(3)(x)  ~   E  -t-^ f(a.) 

i=1T(x-aj} 

J=l 
and 

fU)(x)  "  I      — 2i f(a.) 

1-1  —  (x-&J) 

The  program  then  generates  the  coefficients  of  f(a.)  for  use  by  other  programs 
In  calls  from  XSTP,  the  actual  value  is  desired,  so  the  coefficients  are 
multiplied  by  f (a. )  and  the  value  returned. 
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In  the  case  of  the  Hermite  polynomials,  the  problem  is  to  obtain 
estimates  of  f   (x)  and  f   (x)  when  given  that 

n 

f(x)  =  2  h„(x)  f(a, )  +  h.(x)  f (a.)  (32) 

i=1  i       l     i        i 

For  this  the  derivatives  of  h. (x)  and  h.(x)  are  desired.   The  results  are 
given  by  (assuming  four  points) 

h[6)(x)  =  20(l-a!U.)(x-ai))(il",(x))2+  120(-2£.!(ai))U"(x)  V"U)) 

hfT)(x)  =  lU0(-2£!(a.))(£"'(x))2 
i  'ii 

and 

hj6)(x)  =  20(x-a.)(.r(x))2  +  120^(x)  ^'.'(x) 

hjT)(x)  =  1U0(J?'.'(x))2 

In  the  case  of  the  trimite  polynomials,  the  desired  estimates  are 

f   (x)  and  f    (x).  The  polynomial  f(x)  is  given  by 

n 
f(x)  =  l      (t.(x)  f(a.)  +  t.(x)  f'(a.)  +  t.(x)  f"(a.)         (33) 
i=l 

so  that  the  ninth  and  tenth  derivatives  are  desired  for  t. ,  t. ,  and  t. .   Each 

ii  l 

of  these  has   the   form  V  =  A(£(x))      as   developed  in   section  2.2  where  V  is  one 
of  t. ,   t. ,   or  t..      V^9j(x)    and  V(10'(x)    are  then   given  by 

V(9)(x)   =   (l68o(£*"(x)3)-A  +  15120•Ji"(x)(r"(x))2•A, 

+  36(630U"(x))2  V"(x)    +  U20£'(x)(£"»(x))2)-An 
V(l0)(x)  =16800(2,  (x))3-A»  +  T560QH"(x)U  (x))2-A" 
It  will  be  noted  that  in  all  of  the  cases  above  the  estimates  turned 
out  to  be  functions  of  the  Lagrangian  polynomials.   Therefore,  the  program 
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generates  ^(x).  £  «(x),  £'(a.),  £»(x),  £»(a.),  *"(x)  as  needed  ^  then 

combines  them  to  form  the  desired  estimate. 

The  second  task  of  this  subprogram  is  to  interpolate  to  find 
values  at  new  partition  points.   Since  these  interpolating  polynomials 
are  also  functions  of  the  Lagrangian  polynomials,  much  of  the  code  used 
above  can  also  be  used  here . 

Since  new  partitions  will  only  involve  the  same  or  lower  ordered 
methods.,  the  Lagrangian  method  only  requires  (31)  with  j  =  0.   For  the  Hermite 
polynomials,  the  first  derivative  and  the  functional  values  are  needed. 

The  functional  values  are  obtained  from  (32).   Since  h.  and  h  can  be 

i  i 

written   in  the  form  V  =  A£2(x)   where  A  is    [1-2SL  !  (a.  )  (x-a.  )  ]    and   (x-a.  ), 

respectively,  the  first  derivatives  can  be  written  using 

V'(x)  =  A'£2(x)  +  2M(x)  £»(x) 
Finally,  the  trimite  method  will  require  the  first  and  second  derivatives 
as  well  as  the  functional  values.   The  functional  values  are  obtained 
using  (33).   Using  the  form  V  =  A  3(x)  for  the  fs  (see  page  15  for  values), 
the  first  and  second  derivatives  can  be  obtained  by 

V'(x)  =  A'£3(x)  +  3M2(x)  I  »(x) 
and 

V"(x)  -  A"£3(x)  +  6A'£2(x)  IHx)    +  3A[2£(x)U'(x))2  +  *2(x)  *"(X)] 

Since  calls  to  DERIV  for  purposes  of  interpolating  new  partition 
values  may  require  the  coefficients  in  some  cases  and  values  in  others,  both 
of  these  options  are  made  available. 


52 

3 . h     DIFSUB 

This  program  is  a  modification  of  a  program  by  Gear  [3].  The 
modifications  basically  deal  with  four  areas: 

1.  the  presence  of  the  asymptotic  error  equation, 

2.  the  use  of  DEC0MP  and  S0LVE, 

3.  the  use  of  algebraic  equations  for  the  boundary  conditions, 
h .   the  presence  of  externally  changed  parameters . 

Probably  the  largest  number  of  modifications  deals  with  the  area 
of  the  asymptotic  error  equations.   The  main  item  here  is  that  every  time 
something  is  done  with  the  solution  vector  (Y),  the  same  thing  is  also 
done  with  the  asymptotic  error  system. 

Since  the  solution  to  the  asymptotic  error  system  (HMY)  does  not 
need  to  be  as  accurate  as  Y,  coefficients  were  developed  by  Gear  [k]   so 
that  the  first  two  coefficients  would  be  the  same  in  all  methods,  and  thus, 
allow  the  asymptotic  error  system  to  use  a  lower  ordered  method. 

This  can  be  accomplished  since  the  solution  involves  the 
Newton  iteration 

z  1    j.->\   =  z  /    \  +  &  [-TT—  £]-1  F(z   ,  J  (3*0 

n,(m+l)    n,(m)   —  3a.  —      n,(m)'  VJ  ' 

where  F( z   /  \)  =  hf(u)  -  hu'  in  this  formula, 
n,(.mj       — 

W    L  3a  -J     lZl       M03uJ 
where  I      and  SL      are  the  first  two  coefficients  of  the  method.   Since  W  is 
the  same  for  all  ordered  methods ,  the  two  systems  can  be  solved  using 
different  orders. 
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Another  slight  modification  required  hy  the  asymptotic  error 
system  is  that  in  order  to  start  integrating,  initial  time  derivatives  are 
required.   This  means  that  the  evaluation  of  the  Jacobian  must  be  present 
earlier  in  the  program. 

The  need  for  the  second  modification  is  discussed  in  Gear  [3]. 
He  noted  that  the  break  point  for  using  DEC0MP  and  S0LVE  occurred  when 
about  five  or  more  equations  are  in  the  system.   Since  the  number  of 
equations  solved  here  is  in  general  eight  or  more,  the  need  is  apparent. 
As  an  added  reason,  the  presence  of  the  asymptotic  error  equations  doubles 
the  number  of  calls  to  S0LVE  whereas  the  number  of  calls  to  DEC0MP  will 
remain  about  the  same. 

Two  changes  are  required,  by  the  presence  of  algebraic  equations 
for  the  boundary  conditions.   First,  in  equation  (31)  the  function  F  is 
replaced  by  G(Uj.)  which  is  the  residual  of  the  algebraic  equation.   Second, 
the  calculation  of  those  rows  of  W  involved  with  algebraic  equations  will 
need  to  be  modified.   In  this  case,  the  component  is  given  by 

da  —        0  3a 

where  g  is  the  algebraic  equation.   The  absence  of  I      is  due  to  the  fact 
that  derivatives  of  the  variables  are  not  present  in  G  (see  Gear  [5]). 
In  the  fourth  area,  there  are  two  conditions  that  take  place: 

1.  change  in  the  error  parameter 

2.  change  in  the  partition. 

In  the  first  of  these,  the  only  requirement  is  that  various  error 
condition  values  be  recalculated.   As  a  result,  a  jump  is  made  to  that 
portion  of  the  program. 
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The  second  is  more  involved.   If  the  partition  has  been  changed, 
then  it  will  "be  necessary  to  change  those  variables  dealing  with  the 
number  of  equations  in  the  system,  and  to  re-evaluate  the  Jacobian. 

3.5   DEC0MP  and  S0LVE 

The  purpose  of  these  programs  is  to  solve  Y  =  Ax  for  x  where  A 
is  a  banded  matrix.  The  procedure  is  carried  out  in  two  steps.   (The 
programs  are  modification  of  ones  given  by  Forsythe  and  Moler  [6].)  The 
first  step  is  to  decompose  the  matrix  in  such  a  way  that  the  lower 
triangular  portion  will  represent  the  operations  done  to  the  vector  Y  in 
the  first  half  of  a  Gaussian  elimination  procedure,  and  the  upper  triangular 
portion  will  represent  the  coefficients  of  the  variables  used  in  the  back 
substitution  process.   This  step  is  done  by  DEC0MP  and  is  called  each  time 
the  matrix  A  is  recalculated. 

The  second  step  consists  of  two  parts.   The  first  takes  the 
lower  triangular  portion  and  carries  out  the  operations  of  the  first  half 
of  the  Gaussian  elimination  on  the  vector  Y.  The  second  part  performs 
the  back  substitution. 

In  these  subprograms,  the  matrix  A  is  stored  in  a  double  precision 
array  (UL)  diagonal  by  diagonal,  with  the  main  diagonal  stored  in 

UL(K^+1 ,  I)  where  KT  is  the  number  of  diagonals  and  is  assumed  to  be  odd. 

3.6  DIFFUN 

DIFFUN  calculates  the  time  derivatives  of  the  system  based  on  the 
differential  equation.   It  does  this  by  calculating  the  various  space 
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derivatives  required  and  then  supplies  these  values  to  the  differential 
equation  subprogram  (UFUN).  At  the  same  time,  it  also  calculates  the 
error  estimate  of  the  space  derivatives  by  using  the  stored  coefficients 
of  the  interpolating  polynomials  for  such  estimates.   These  results  are 
then  added  to  the  derivatives  and  given  again  to  UFUN.  The  difference 
between  the  results  of  the  two  calls  to  UFUN  then  provides  an  estimate 
of  dx(xjL,t). 

A  second  use  of  DIFFUN  is  to  help  determine  the  Jacobian  used  in 
DIFSUB  in  the  predictor-corrector  method  and  in  the  differential  equation 
for  the  asymptotic  error.   In  this  environment,  the  error  estimates  are 
not  used  so  a  variable  is  set  which  causes  this  calculation  to  be  bypassed, 

3.7  DIVERR 

The  asymptotic  error  equation  was  given  by  (section  2.2) 
e»(t)  =  f  (y(t))-e(t)  -d-d 

y  t   x 

This  program  calculates  the  time  derivative  e'(t)  for  each  variable  using 
fy(y(t))  which  is  stored  in  the  array  PEQ,  d  which  is  given  by 

ERR0R(IA) 
PERTST   '  and  dx  which  is  st°red  in  ERR(IA). 
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k .      RESULTS 

The  system  was  used  successfully  to  examine  three  separate 
examples.   In  each,  case,  the  problem  was  solved  using  the  system  and  an 
alternative  numerical  technique.   The  three  problems  were: 

1.  u  =  u   ,  u(0,t)  =  u(l,t)  =  0, 

t     XX 

u(x,0)  =  sin(ut)  (35) 

2.  ux  =  u   +  u  +  u,  u(0,t)  =  u(l,t)  =  0, 

t     XX     X 

u(x,0)  =  e"X  sin(.TTt)  (36) 

3.  u  =  u  ,  u(0,t)  -  u  (0,t)  =  0 

t      XX  X 


u(l,t)  +  u  (l,t)  =  0,  u(x,0)  =  sin  (ttx)  (37) 


Example  1 

In  this  example,  the  solution  is  quite  simple  to  calculate  and  is 


given  by 


-TT2t 
u(x,t)  =  e~    sin(iTx)  (38) 


The  independent  solution  uses  a  method  developed  by  Crank  and  Nicholson  [Tl 
The  programming  of  the  method  is  based  on  information  in  Richtmyer  and 
Morton  [8].   The  program  is  listed  in  Appendix  B. 

In  solving  the  problem  with  the  latter  method,  an  answer  was 
obtained  using  a  step  size  of  Ax  =  g  and  At  =  Ax.   The  process  was  then 
repeated  with  step  sizes  of  -7-,  r^",  ZJ7»-"  until  the  error  was  less  than 
.001.   Results  are  in  Table  h. 
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In  using  the  system,  answers  were  obtained  with  EPS  equal  to 
.1   .01,  and  .001.  The  errors  in  the  first  two  cases  are  approximately 
the  same  since  the  same  initial  partition  was  used  in  both.   The  results 
are  given  in  Tables  5,6,  and  7.   In  these  tables,  the  column  headed 
"Number  of  Equations  Evaluated"  represents  the  number  of  calls  to  UFUN 
times  the  number  of  equations  per  partition  point. 

The  solution  for  EPS  =  .001  when  using  the  system  took 
approximately  26  seconds.   If  it  had  been  known  that  Ax  =  ^  would  have 
produced  the  desired  solution,  about  seven  seconds  would  have  been 
required.   It  would  be  normal  in  developing  the  latter  solution  to  compare 

the  results  obtained  with  two  different  step  sizes  to  estimate  the  error. 

1      1 
If  this  were  done  using  Ax  =  ^  and  gjf,  the  difference  would  have  been 

2  2 

about  .002.   Since  the  truncation  error  is0[(At)  ]  +0[(Ax)  ]  [8],  it 

would  be  expected  that  the  error  in  the  solution  using  Ax  =  jpj-  would  be 

near  .0005  which  would  be  within  the  .001  desired.   This  process  would  have 

required  about  eight  seconds. 

If  the  difference  between  the  two  solutions  themselves  must  be 
less  than  .001,  then  step  sizes  of  -^  and  j^q   would  have  been  required. 
Using  these  two  step  sizes,  the  result  would  have  required  approximately 
32  seconds. 

The  first  method  for  obtaining  the  solution  would  require  -  as 
much  time  as  that  obtained  from  the  system.   The  second  method  would  have 
required  more  time  than  the  system  to  obtain  the  solution. 
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Example  2 

Although,  the  equation  (36)  can  be  put  into  the  form  of  equation 

(35)  "by  the  substitution  u  =  e~  v,  it  was  felt  an  example  involving  other 

terms  would  he  desirable.  Because  of  the  transformation  and  (38),  the 

solution  to  (36)  is  seen  to  be 

2, 
u(x,t)  =  e   e     sin(-rrx) 

The  program  used  to  obtain  the  independent  solution  is  listed  in 
Appendix  C.   The  results  of  this  program  are  given  in  Table  8.  The  results 
obtained  from  the  system  for  EPS  =  .1,  .01,  and  .001  are  in  Tables  9,  10, 
and  11,  respectively. 

For  the  independent  solution,  a  step  size  Ax  =  ^TZ*»  At  =  Ax  was 
required  before  the  solution  was  in  error  by  no  more  than  .001.  However, 
the  largest  difference  between  solutions  with  Ax  =  7^,  Ax  =  r^g 
was  .0015  which  would  indicate  that  Ax  =  tj-^q-  could  be 

used.   If  Ax  =  -77-  and  Ax  =  — ^  were  used  to  determine  the  solutions,  31 
seconds  would  be  used.   If  Ax  =  rrr,  and  Ax  =  ^were  used,  then  12U 
seconds  would  be  required.   The  solution  using  EPS  =  .001  with  the  system 
also  requires  about  31  seconds. 


Example  3 

The  system  can  also  be  used  to  solve  problems  where  the  boundary 
conditions  are  of  a  more  general  form.   In  this  example,  the  boundary 
conditions  are  u(0,t)  -  u  (0,t)  =  0  and  u(l,t)  +  u  (l,t)  =  0.   The  solutions 
are  expected  to  take  longer  because  the  errors  introduced  at  the  boundary  due 
to  extrapolation  will  require  smaller  step  sizes  in  DIFSUB.  Because  of  this 
only  .1  and  .01  will  be  considered  for  values  of  EPS.  These  results  are  giver 
in  Tables  12  and  13,  respectively. 
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The  independent  solution  was  again  obtained  by  using  the 
Crank-Nicolson  method.   The  program  used  to  solve  this  system  of  equations 
was  based  on  information  from  Ozisik  [9].   The  program  listing  is  given 
in  Appendix  D.   The  results  are  given  in  Table  lk . 

To  obtain  a  solution  using  the  CrankJJicholson  method  requires 
approximately  six  seconds.  A  similar  solution  requires  about  3U  seconds 
using  the  system. 

In  examples  2  and  3  above,  it  was  found  that  results  were 
easier  to  obtain  if  the  auxiliary  algebraic  equations  at  the  boundaries 
were  bypassed.  This  is  done  by  setting  up  equations  of  the  form 

\   =B(u.,  v.,  x.,  t.) 

r=  v.  -B(u.,  v.,  x.,  t.) 
where  B  represents.,  an  interpolation  polynomial.  ' 

In  order  that  the  Jacobian  will  remain  invertible ,  a  1  is 
inserted  in  the  main  diagonal  of  the  equations  bypassed.   The  dependence 
that  interior  points  would  have  had  on  these  boundary  values,  as  present 
in  the  Jacobian,  is  now  transferred  to  the  interior  points  used  to  obtain 
that  value. 

Although  the  error  estimate  used  for  the  spatial  truncation  error 
appears  to  be  accurate  enough  for  the  previous  examples,  the  presence  of  a 
more  accurate  one  may  be  necessary  in  more  complicated  problems.   The 
system  was  used  to  try  to  integrate  a  problem  where  the  initial  value 
involved  a  peak  in  the  center  of  the  interval .  Although  the  system  was  able 
to  produce  values  of  adequate  precision,  which  would  appear  to  indicate  an 
analog  of  sufficient  accuracy,  the  error  estimates  were  large  enough  that 
the  integration  would  have  been  rejected.   (See  Table  15) 
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72 
A  similar  situation  occurred  with  the  following  example: 

Ut  =  87  (l+u)  fp  u^°»t)  =  u^1»t)  =  °  u(x,0)  =  sin3(irx) 
x  e[0,l] 
For  use  with  the  system,  the  equation  was  rewritten 

Ut  "  (3^   +  (1+U)  72 

dX 

An  independent  solution  was  obtained  by  using  a  predictor- 
corrector  modification  of  the  Crank-Nicholson  method  that  has  been  developed 
by  Douglas  and  Jones  [10],   The  solution  was  also  obtained  for  a  few  values 
of  t  using  the  system  with  the  error  control  portions  bypassed.   The  results 
agreed  with  those  obtained  with  the  alternative  method.  However,  the  error 
estimates  for  the  space  analogs  would  have  caused  the  integration  to  cease. 
(See  Table  l6) 

The  results  obtained  from  the  first  three  examples  appear  to 
indicate  the  feasibility  of  using  the  method  of  lineswith  the  subprogram 
DIFSUB  to  produce  an  automatic  program  when  the  asymptotic  error  equations 
are  integrated  along  with  the  original  system.   The  two  additional  examples 
seem  to  provide  the  possibility  of  integrating  other  types  of  problems 
including  those  of  a  nonlinear  nature. 
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APPENDIX  A 
THE  SYSTEM  OF  PROGRAMS 
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SUBROUTINE  P ARDIF (R INT, NCK, DI VER, NPS AVE , PS AVE , EPS ,SMAX, SMI N, SX,  PDF  1000 

1SAVE,NEQ,MAXAR,MAXN, IDE R, MS TRT , INP,NUM, MOLD , Y, X,H,HMIN,HMAX , YMAX ,  POFlOiO 

2ERR0R,KFLAG,JSTART,MAXDER,K7,HMY, EROR,T, I  NT ,HSAV , PEQ, ERORA I  PDF  10  20 

C  PDF1030 

C  THE    PURPOSE   OF    THIS    PROGRAM    IS    TO    INTEGRATE    ONE    STEP,       TO    DO  PDF1040 

C  THIS    IT    CREATES    A    PARTITION,     INTEGRATES,    AND    ATTEMPTS    TO    CONTROL  PDF1050 

C  THE    ERRORS  PDF1060 

C  THE    VARIABLES    ARE  PDF1O70 

C  RINT  BOUNDARY    POINTS  PDF1080 

C  NCK  USED    IN    XSTP  PDF1090 

C  DIVER  CONTAINS    ERROR    EXTIMATE    COEFFICIENTS  PDF1100 

C  NPSAVE  USED    IN    XSTP  PDF1110 

C  PSAVE  CONTAINS    DECOMPOSED    JACOBIAN  PDF1120 

C  EPS  ERROR    PARAMETER    -    USER  PDFU30I 

C  SMAX,  PDF1140I 

C  SMIN  LIMITS    ON    SPACE    STEP    SUE  PDF115>0 

C  SAVE  STORES       Y  PDFU60I 

C  MSTRT,  PDFUTOi 

C  MAXAR  LIMITS    ON    ANALOG    TYPE  PDF1180S 

C  1    -    LAGRANGIAN  PDF119Q) 

C  2    -    HERMITE  PDF1200? 

C  3    -    TRIMITE  PDF1210( 

C  IDER  DERIVATIVES    USED  PDF1220( 

C  Y  ANSWER    VECTOR  PDF123QI 

C  X  PARTITION  PDF1240I 

C  HMIN,  PDF125Q* 

C  HMAX  LIMITS    ON    TIME    STEP    SIZE  PDFl260t 

C  YMAX  USED    FOR    RELATIVE    ERROR    WHEN    Y    IS    GREATER    THAN    1  PDF1270C 

C  ERROR  TIME    ERROR  PDF1280( 

C  KFLAG  RETURN      CODE  PDF1290( 

C  HMY  ASYMPTOTIC    ERROR    VALUE  PDF1300C 

C  EROR  SPACE    ERROR  PDFi310< 

C  HSAV  SAVES    HMY  PDF1320( 

C  PEQ  JACOBIAN  PDF1330C 

C  ERORA  ESTIMATED    ERROR    INTRODUCED    IN    THIS    STEP  PDF1340( 

C  PDF1350C 

IMPLICIT    REAL    *    8(A-H,P-Z)  PDF1360C 

DIMENSION    ERORAdl  PDF137QC 

DIMENSION    KL8(3I,RINT(2),NCK(1I ,D I VER (9, 1 ) , NPSAVE (1  I, PSAVE (1) ,SX<  3PDF1330r 

1),NEQ(3»,IDER<4,3), INT( 2, 1 ) , X{ i I , Y( 5 , 1) ,YT(3I  ,YTA ( 3 , 5 ) , XT ( 5  I  ,  PDF139K 

2ERROR(l>,EROR(l ) ,  YMAX  { 1  ) ,  HMY  (  2, 1  > ,GAM( 3 ) , SAVE (7 , 1  I ,HSAV( 3 ,1 ) ,  PDF1400C 

3PEQ(1 »,MAXM(3}  PDF1410C 

COMMON    /DVR/    DVY( 41  ,COEF( 12     , 3 » , ALPHAi 2, 4,3 » , DI V( 16 , 5  I  PDF1420C 

COMMON    /CNT/    ICSTPS ,NUMDE ,MUMSO  PDF1430C 

DATA       PVAL/.5/  PDF14'?0C 

DATA    RVAL,SRVAL/. ID0,.01D0/  PDFl4iJi 

DATA    KL8    /9,19,29/  PDF1460C 

PDF147GC 
AND   CHECKS    FOR    CHANGE    CONDI TIONSPDF148D0 

POF149O0 
PDF15000 
PDF15100 
PDF15200 
PDF15300 


c 
c 
c 

2 

THIS  SECTION  INITIALIZES  COUNTERS 

IFUSTART  .NE.  0  1   GO  TO  10 
MAXSV  =  MAXAR 
ISPSCT  =  0 
TPS  =  .4  *  EPS 
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SPS  »  .1  *  EPS 


DO  15  IA*1,3 


INP  *  1 
SEPS  =  EPS 
JSTART  »  0 
KTPSCT  =  0 
ICPART  =  0 
KSPSCT  =  0 
ITUPS  =  1 


10   IF(SEPS  .EQ.  EPS)   GO  TO  23 


TPS  =  .4  *  EPS 


INP=2 
J=NEQ(MO) 


PDF15400 


VEPS  *  DLOGCEPSI  POF15500 

Nl«7*MAXN*l  POFISftOO 


PDF15800 


15   MAXMCIA)  *  MAXN/IA  -2  pnF  „;; 

NUMA«MAXN  PDF15900 


P0F16000 
PDF16100 
PDF16200 
POF1630O 
PDF16400 
P0F16500 
P0F16600 
PDF16700 


J^V  X  P0F16800 

iC     U  PDF16900 

NCNGE  '    °  pnPWAon 

ITPSC   -  0  0F17000 

ICSTPS  *  0 
NUMDE  «  0 
NUMSO  «  0 
T  »  0.000 
MAXER  =  MAXSV 
GO  TO  200 


PDF17100 
PDF17200 
PDF17300 
P0F17400 
PDF17500 
PDF17600 

20   IFCNCNGE  .EQ.  0»   GO  TO  255  PDF17800 

NCNGE  "  °  P0F17900 

GO  TO  12 


PDF18000 
PDF18100 


^'S  =  ?L°G1^S>  PDF18200 


PDF18300 


SPS=.1*EPS  pnF1    " 

12   IFilCSTPS.LT.  101   GO  TO  5  PDF18500 

200   CALL  XSTPCRINT.NCK, DIVER, NPSAVE I Nl I , NPS AVE ( N2 ) , SPS, SMAX,SMIN,SX  PDF18600 

1»KFLAG,NEQ,MAXER,MAXM,IDER,MSTRT,PSAVE,INP,NUMA,M0LD,X,Y,M0)  PDF18  700 

IFCKFLAG  .LT.  C   .AND.  MAXSV  .EQ.  1>  RETJRN  PDF18flOO 

IFCKFLAG  .LT.  0)     GO  TO  299  PDF18900 

IFCICK  .EQ.  1   .AND,  MO*NEQ(MOI     . GT.  NUMI   GO  TO  208  PDF19000 

MAXER  =  M0  PDF19100 

THIS  SECTION  PUTS  THE  NEW  PARTITION  INTO  USABLE  FORM  AND  PDF19300 

OBTAINS  THE  FUNCTIONAL  VALUES  FOR  THIS  NEW  PARTITION.  PDF19V30 

PDF19500 

PDF19620 

,tk,T.,»  PDF19700 

52-Soi?  li7?I?  PDF19800 

JA-MO*( J+l»*l  PDF19900 
SAVE(1,JA                       I-RINTC2I 
JA=M0-1 

IFCJA.LE  .01  JA=3 


PDF20000 
PDF2C100 

IFCJSTART.EQ.OI  INT(2,l)«l  PDF23?0ft 

DO  222  IA=1,J  PDF20500 

N3=  M0*3*(IA-1)4N2-1  PDFPD600 

IB=NPSAVECN3I  PDF20600 

id  i»ri«vi:i(N3i  PDF2070n 

INT(1,IA*U  =  DIVERCM0,IB)  PDF2C800 
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222 


C 
C 

c 


250 

208 


300 


IAA»l*MO*(IAI 

SAVEl  1,  IAA) -DIVER  I  MOO,  IB  I 

N3»N3-MOOA 

IF( IA.EQ.l»NPSAVE(N3»»l 

N3»N30 

IFUA.EQ.l)NPSAVElN3)»l 

IFl IA  .EQ.  2)  NPSAVEIN3I  »  2 

IFIIA  .GT.  2   .AND.   IA  .LE.  Jl  NPSAVECN3) 

IF(  IA.EQ.  J.  AND.MO.EQ.il  NPS  AVEl  N3  )»J-2 

IFUSTART.EQ.OJ     INT (2, I A+l ) »NPSAVE(N3> 

CONTINUE 

IFUSTART.NE.OI    NPS  AVE!  N30  l»NPSAVE(N3) 

NUMA=NEQ(MOJ 

IFl JSTART.GE.l)    GO   TO   300 

NEW  VALUES  FOR  INITIAL  PARTITION. 

IQQQOO 

J2»JO 

DO  250  IA*1,J2 

IAA1  «  MO*( IA-1  I 

IAA  =  IAA1  +1 

XI IA»=SAVEI1,IAA» 

CALL  INITIYT,M0,XUAI> 

00  250  IB=l,MO 

IAA  »  IAA1  ♦  IB 

HMY(1,IAA)  »  0.0 

Yll,IAA)=YT(IBI 

GO  TO  258 

MO  =  MOLD 

J  *  NUM 

NEOIMOI  =  J 

GO  TO  258 


310 

10101 

350 
C 
C 


IA  -I 


/MOLD 


THIS  SECTION  FIRST  CHECKS  TO  SEE  IF  THE  POINT 
PARTITION.  IF  IT  WAS,  THE  VALUES  ARE  REUSEO. 
VALUES  ARE  OBTAINED  BY  INTERPOLATION. 


IW=JSTART+1 

NUM=NUM/MOLO 

J2=J  +2 

DO   3  90  IB=1,J2 

IA  =  1 

IAA=l+M0*(IB-ll 

XI=SAVEll,IAA» 

IFIXI.LE.X1  IA  I  I 

IA=IA+1 

IFl IA.LE.NUM+2) 

KFLAG=-30 

RETURN 

KFLAG=-41 

RETURN 

IFl  XI.EQ.XC  IAII 


WAS  IN  THE  EXISTING 
IF  NOT,  THE  NEW 


GO  TO  350 


GO  TO  310 


GO  TO  370 


INTERPOLATION  SECTION 
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KG»0  PDF26400 

351  IF<  IA.LT.NUM*2HC  =  INT(2,IA)  PDF26500 
IF(  IA.  EQ.NUM+2)  IC=NUM-1  P0F2660O 
IW1  «  I  PDF26700 

352  DO  355  IE-1,4  PDF26800 
ID«IC-1*IE  PDF2690O 

355   XT(IEI«X(IOI  PDF27000 

DOAL35fir=XlT:MOA'IDER,l,MO,,MO,1'l'KG'MC,LD»XI'GAM'XHI  pSfItISS 

IAB=IE4M0*( IB-ll  PDF2730O 

IF(IW1  .EQ.  II  SAVE(IW2,IAB>  »  0.0  oncflc?? 

IFdWl  .EO.  1)  HSAV(l.IAB)  »  0.0  n°      ? 

HSAV<2,IAB)  =  0.0  POF27600 

DO  358  ID  =  lUlftM  PDF27700 

SAVE! ID*ltIABI»0.  P0F2r800 

DO  358  IF=l,MOLD  PDF27900 

DO  358  IG=l,4  PDF28000 

IH-IC-1+IG  PDF28100 

IAA=IG+4*(IF-1»  PDF28200 

IAC=IF+M0LD*(IH-1»  PDF28300 

358      >D  :^ 

IFIIC  .ST.  I   .OR.   IC  .LT.  NUM-1,   GO  TO  390  £™0° 

IF(IC  .EQ.  1»  IC  =  2  PDF29000 

IFIIC  .EQ.  NUM-ll  IC  =  IC  -1  PDF29100 

GO  TO  352  PDF29200 

:  PDF29300 

USING  THE  OLD  PARTITION  VALUES  an™™ 

370   DO  371  IJ  =  ltMO  PDF29600 

IAA  =  I J  +  M0*( IB-1  I  PDF29700 

IA3=IJ+M0LD*(IA-1»  PDF29800 

SAVE(IW*2,IAAJ=YMAX(IAB)  PDF29900 

HSAV<1,IAA|  =  HMYU.IABI  PDF30000 

HSAV(2,IAAJ  =  HMY(2,IA8)  PDF30100 

DO  371  10=1, IW  PDF30200 

SAVE<  ID-H,IAA»  =  Y<  ID.IAB)  PDF30300 

CONTINUE  PDF30400 

PDF30500 

UPDATES    THE    VECTORS    FOR    THE    NEW    PARTITION  PDFSO^OO 

J2=J*2  PDF3CB0O 

DO      410    IA=1,J2  PDF30900 

N3=N2-l*JA+3*( IA-ll  PDF31000 

INT(2,IAI=NPSAVE(N3)  P0F31100 

IAB=1*M0MIA-1>  PDF31200 

X(IA»=SAVE{1,IAB)  PDF31300 

DO    410    IC=1,M0  PDF31400 

IAA  =  IOM0*(  IA-1  I  PDF31500 

YMAX(IAA»=SAVE(  IW2.IAAJ  oS«fJ2! 

DO    41C     IBM.IW  PDF31700 

PDF31800 
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HMY(1,IAAI     »    HSAV(1,IAAI  PDF3190( 

HMY(2,IAAI    »    HSAV(2,IAA)  P0F3200C 

410       Y(IB,IAAI=SAVE(  IB+1  ,IAA)  PDF3210( 

DO    430    IA    *    1,4  PDF3220C 

430      XT(IA)    »    XllA  +  ll  POF3230C 

XI    «    X(l>  PDF3240C 

CALL    DERIV(XT,YTA,IDER(l,MOI,l,l,l,0,l,XI,GAM,XHI  POF3250C 

DO    435    IA    =    ltMO  PDF3260C 

DO    435    IB    *    2,IW  PDF3270C 

YUBfIA)    =0.0  PDF3283C 

IF(  IB    .EQ.    21    HMY(2,IA)    »    0.0  PDF3290C 

DO    435    IF    *    lf4  PDF3300C 

IE    =    IA    ♦    MO    *( IF-1  I  PDF3310C 

Y(IB,IAI    =    Y(IB,IA»    +    COEF< IF,1)*Y(IB,IEI  PDF3320C 

IF(     IB    .EQ.    21    HMY(2,IA»    ~    HMY(2,IA)    ♦    COEF  ( I  F,  1  l*HMY  (  2,  IE  I  PDF3330C 

435      CONTINUE  PDF33400 

DO    440    IA    *    1,4  PDF3350G 

440       XT(IA)    =    X(J2    -     IA>  PDF3360C 

XI    =    XU2I  PDF3370C 

CALL    DERIV<XT,YTA,IDER(1,M0I, 1,1, 1,0,1, XI, GAM, XHI  PDF33800 

DO    445    IA    =    1,M0  PDF33900 

DO    445    IB    *    2,IW  PDF34000 

IG    =    J2    -    MO    ♦    IA  PDF34100 

Y(IB,IGI    =    0.0  PDF34200 

IF( IB    .EQ.    2)    HMY(2,IG>    »    0.0  PDF3430Q 

DO    445    IF    =    1,4  PDF34400 

IE    ■    IA    ♦    MO    *    (J2    -    IF    -II  PDF34500 

YU8,IG>    =    Y(IB,IG)    ♦   C0EF(IF,1»    *    Y(IB,IEI  PDF34600 

IF(IB    .EQ.    21    HMY(2,IG)    =    HMY<2,IGI*    COEF( IF, 1 )*HMY(2 , I  El  PDF34700 

445      CONTINUE  PDF34800 

IDOUB=IH  PDF34900 

IQQQC=1  PDF35000 

C  PDF35100 

C      CALCULATES  THE  COEFFICIENTS  FOR  THE  ERROR  ESTIMATING  INTERPOLATINGPDF35200 

C      POLYNOMIAL  PDF35300 

C  PDF35400 

258   DO  262  IA=1,J  PDF35500 

ID=MO*IA«-l  PDF35600, 

XI=X(IA*1I  PDF35700 

DO  260  IB=1,5  PDF35800 

IC=INT(2,IA  +  1I-1+IB  PDF35900 

XT(IB>=X(ICI  PDF36000 

260      CONTINUE  PDF36100 

IAA=INT(1,IA*1  I*  I  A  PDF36200 

XH=XI~X(IAA<-1  I  PDF36300 

CALL    DERIV(XT,YTA,IDER(1,M0|, MO, 1,0,1, MOLD, XI, GAM, XHI  PDF364Q0 

DO    262    IB=1,M0  PDF36500 

IAA1    =    9*( IB-ll  PDF36600 

DO    262    IC=1,9  PDF36700 

IAA    =    IAA1    ♦    IC  PDF36800 

262      DIVER( IAA,IDI=COEF( IC,IBI  PDF36900 

K7=KL8(M0I  PDF37000 

N=(NEQ(M0»*2I*M0  PDF37100 

C  PDF37200 

C               INTEGRATES    AND    DETERMINES    THE    MAGNITUDE    OF    THE    ERROR  POF37300 
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255 


270 


275 


288 


291 


CALLDIFSUB(N,T,Y, SAVE, H,HMIN,HMAX, TPS, YMAX, ERROR, KFLAG.JSTART. 

2hSq!S 
ICSTPS    *    ICSTPS    ♦    1 
ICK    »    0 

IQQQC=0 

NUM=NEQ(M0)*M0 

MOLD* MO 

IWJSTART  +  1 

IF(KFLAG    .LE.    -3)       GO   TO    290 

IF(KFLAG.LE.O)    RETURN 

NUM1  =  (MUM*M0 

NUM2=M0*1 

RMAX    =   0.0 

DMAX      =    0.0 

00    27C    IA=NUM2,NUM1 

DBStA    =    DABS(E    RORA(tAi) 

IF(DBSEA.GT.RMAX)    RMAX    =    DBSEA 

IF{ ICPART    .LT.    31       GO    TO    270 

DMAXT    =    DABSIEROR(IA) I 

IF(DMAXT      .GT.       DMAX)       DMAX    =    DMAXT 

CONTINUE 

CHECKS  ON  ACTIONS  BASED  ON  ERRORS  AND  COUNTERS 
IF(RMAX.GE.EPS)  GO  TO  290 
ERROR  IS  SMALL  ENOUGH 

SEPS    =    EPS 

IF(RMAX    .GE.    RVAL*EPS    . OR. ( ICSTPS-I RECSP >    .LT.     10)    RETURN 

ICPART      =    ICPART    ♦    I 

IRECSP    =    ICSTPS 

IF( ICPART    .LE.    31       GO    TO    275 

ICPART    =    0 

IFIDMAX    .LT.     SRVAL*EPS)       GO    TO   292 

IF    (ITUPS    .GT.    3)       GO   TO    288 

CHANGES    TIME    PARAMETER 

VTPS  =  DLOG  (  TPS  ) 

VTPS  =  (  VTPS  ♦  VEPS  I  /  2. 

TTPS  =  DEXP  (VTPS) 

TPS    =    DMIN1        (TTPS, TPS*. 5 DD* ( EPS-RMAX > ) 

ITUPS    =    ITUPS    ♦    1 

RETURN 

IF(RMAX    .GT.     SRVAL*EPS)       RETURN 

CHANGES    SPACE    PARAMETER 

IF(SX(MO)*2.0DC    .GE.    SMAX       .AND.       MO    .EQ.     MSTRT)       RETURN 
ICK     =     1 

TSPS    =    2.    *    SPS 

SPS   =  DMINKTSPS,  .5  *  EPS) 


PDF37400 

PDF37500 

PDF37600 

PDF37700 

PDF37800 

PDF37900 

PDF38000 

PDF38100 

PDF38200 

PDF38300 

PDF38400 

PDF38500 

PDF38600 

P0F38700 

PDF38800 

PDF38900 

PDF39000 

PDF39100 

P0F39200 

PDF39300 

PDF39400 

PDF39500 

PDF39600 

PDF39700 

PDF39800 

PDF39900 

PDF40000 

PDF40100 

PDF40200 

PDF40300 

PDFA0400 

PDF40500 

PDF40600 

PDF40730 

PDF40800 

P0F40900 

PDF41000 

P0F41100 

P0F41200 

PDF41300 

PDF41400 

PDF41500 

PDF41600 

PDF4170Q 

PDF41800 

PDF41900 

PDF42000 

PDF42100 

PDF42200 

PDF42300 

PDF42400 

PDF42500 

PDF42600 

PDF42  700 

PDF428Q0 
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292 


290 


296 


298 


299 


29  3 


295 


294 


TPS      »    .6    *    EPS      -SPS 

IF(  SX(M0I*2.000    .GE.    SMAX      .AND, 

ITUPS    »    I 

ICPART    *    0 

NCNGE    «    I 

RETURN 

ERROR    TOO   LARGE 

ITUPS    =    I 

ICPART    «    0 

ICSTPS    =    ICSTPS    -    1 

IRECSP    =    ICSTPS 

IF(  ICSTPS    .LT.    101    GO    TO    298 

IF(KFLAG    .LE.    -3»       GO    TO    293 

IF(TPS    .LT.    .00001    .OR.    KTPSCT   < 

GO    TO    294 

KSPSCT   =    KSPSCT    ♦    I 

KTPSCT    »1 
IF(KSPSCT    .LE.    101    GO   TO   293 
KFLAG    =    -42 
RETURN 


IF    THE    NUMBER   OF    STEPS    IS    SMALL    ENOUGH 

PARAMETERS 

SPS    =    .1    *    SPS 

TPS    =    .1    *    TPS 

NCNGE    =    1 

ISPSCT    =    ISPSCT    ♦    1 

IF(  ISPSCT    .LT.4    I    GO   TO    10 

KFLAG    =   -35 

IF(MAXSV    .EQ.    1  J      RETURN 

MAXSV    =    1 

IF( ICSTPS    .GT.    10    I       RETURN 

GO    TO    2 

DECREASES  SPACE  PARAMETER 


MO  .EQ.  MSTRTI   RETURN 


GE.  5  I  GO  TO  296 


-  RESTARTS  WITH  SMALLER 


GO  TO  295 


IF(TPS  .LT.  .00001) 

TPS  =  .000001 

J  START  =  -I 

GO  TO  255 

SPS  =  PVAL  *  SPS 

GO  TO  400 


DECREASES  TIME  PARAMETER 

VTPS  =  DLOG(TPS) 

VTPS  =  (VTPS  ♦  DL0G(i.0-6M/2.00 

KTPSCT   =  KTPSCT  ♦  1 

TPS  =DEXP  {  VTPS) 

ITPSC  =  i 

JSTART  =  -1 

GO  TO  255 


PDF42900 
PDF43000 
PDF43100 
PDF43200 
PDF43300 
PDFA3400 
PDF43500 
PDF43600 
PDF43700 
PDF43800 
PDF43900 
PDF44000 
PDF44100 
PDF44200 
PDF44300 
PDF44400 
PDF44500 
PDF44600 
P0F447C0 
PDF44800 
PDF44900 
PDF45000 
PDF45100 
PDF45200 
PDF45300 
PDF45400 
PDF45500 
PDF45600 
PDF45700 
PDF45800 
PDF45900 
PDF46000 
PDF46100 
PDF46200 
PDF46300 
PDF46400 
PDF46500 
PDF46600 
PDF4670Q 
PDF46300 
PDF46900 
PDF47000 
PDF47100 
PDF47200 
PDF473G0 
POF47400 
PDF47500 
PDF47600 
PDF47700 
PDF47800 
PDF47900 
PDF48000 
PDF48100 
PDF48200 
PDF483D0 
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c 

OBTAINS    THE    OLD    PARTITION 

c 
c 

400 

CALCULATED 

CONTINUE 

DO    405    KF1»1,N 

HMY(l,KFl»    *    HSAV<l,KFl) 

HMY(2,KF1I    «    HSAV(2,KF1) 

DO    405    KF2=1, IW 

405 

Y(KF2,KFI»=SAVE(KF2,KF1» 

GO    TO    200 

END 

SO  THAT  THE  NEW  PARTITION  MAY  BE 


PDF48400 
PDF48500 
PDF48600 
PDF48700 
PDF48800 
PDF48900 
PDF49000 
PDF49100 
PDF49200 
PDF49300 
PDF49400 


SUBROUTINE  XSTP (R INT , NCK, DPOS ,NPT .NPTER .ERQR. HMAX  hmtw  cv  *    Men 


THE  PU 

PARTIT 

THE  VARIABLES 

RINT 

NCK 

DPOS 

NPT 

NPTER 

EROR 

HMAX 

HMIN 

SX 

K 

MAXO 

MAXM 

IDER 

MSTRT 

YTEMP 

INP 

NUM 

MOLD 

X 

Y 

MONW 

AF 

AD 

AVA 

AA 
AE 
YTfXT 

L 


RPOSE  OF  THIS 
ION  OR  MODIFY 


ITS 


THROUGH 


SUBPROGRAM  IS  TO  PROVIDE  AN  INITIAL 
AN  EXISTING  ONE 
IN  THIS  SUBPROGRAM  ARE 

BOUNDARY  POINTS 

AN  ARRAY   +1  INDICATES  POSITION  SATISFIES 

ERROR    CONDITION;    -1       INDICATES    IT    DOES    NOT 

CONTAINS  PARTITION  POINTS  ALONG  WITH 

ANALOG  NE IGHBORS 

POINTER  ARRAY  FOR  OLD  £  NEW  PASSES 

THE  INTERVAL 

OUTPUT  POINTER  FOR  ARRAY 

SPATIAL  ERROR  PARAMETER 

LARGEST  ALLOWABLE  STEP  SIZE 

SMALLEST  ALLOWABLE  STEP  SIZE 

CONTAINS  THE  SMALLEST  STEPS  USED 

RETURN  CODE 

MAXIMUM  ORDER  ALLOWED 

MAXIMUM  NUM8ER  OF  POINTS  IN  EACH  ORDER 

CONTAINS  COOE  FUR  DERIVATIVES  USED 

LOWEST  ORDER  ALLOWED 

TEMPORARY  STORAGE  OF  PARTITION  FUNCTION 

1  -  INITIAL  PARTITION;  2  -  MODIFY 

NUMBER  OF  POINTS  IN  OLD  PARTITION 

OLD  ORDER 

OLD  PARTITION 

OLD  FUNCTION  VALUES 

NEW  ORDER 

TEMPORARY  STORAGE  FOR  TIME  DERIVATIVE 

TEMPORARY  STORAGE  FOR  ERROR  ESTIMATE 

RELATIVE  POSITION  OF  POINTS  USED  IN 

PARTITION  CALCULATION 

TEMPORARILY  STORES  SPACE  DERIVATIVE 

TEMPORARILY  STORES  PERTURBATION  OF  AA 

STORES  VALUES  FOR  DERIV 

RECORDS  POSITIONS  FOR  POINTS  IN  ERROR  EXTIMATE 


VALUES 


IN  AF 
INITIAL 


STP10000 

STP10100 

STP10200 

STP10300 

STP10400 

STP10500 

STP10600 

STP10700 

STP10800 

STP10900 

STP110C0 

STP11100 

STP11200 

STPH300 

STP11400 

STP11500 

STP11600 

STP11700 

STP11800 

STP11900 

STP12000 

STP12100 

STP12200 

STP1230O 

STP12400 

STP12500 

STP12600 

STP12700 

STP12800 

STP12900 

STP13000 

STP13100 

STP13200 

STP13300 

STP13400 

STP13500 

STP13600 

STP13700 

STP13800 
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FOR  DVY, ALPHA, DIV,COEF, GAM        SEE  DERIV 


10 


20 


30 


IMPLICIT  REAL  *  8  (A-H,CW) 

DIMENSION  RINTI2) ,NCK(1I , DPOS (3 , 3 , 1 ) ,NPT < 2, 1 > ,NPTER (3 ,1 ) , A( 3, 2 , 5) 
1  AFC  3)  ,AD(3) ♦ GAM  1 3  ) , AVA ( 5 » , A  A ( 4) , AE( 4) ,HQ( 21 »HSQ< 2 ) , IDER <4 , 3 ) , 
2YT(3,5),XT(5),YTEMP(3,1),L<5) 

OIMENSION  SX<3),NEQ(3) 

DIMENSION    ALFA<2) , X ( 1 ) , Y< 5, MOLD, I » ,MAXM( 31 

DIMENSION    I  WD  1(3,4) ,R WD( 3 ,2 ) , WVAR ( 3 ) 

COMMON    /CNT/    ISTEP, NUMDE , NUMSO 

COMMON    /DVR/    DVY( 4 >  ,COEF ( 4, 3 , 31 , ALPHA ( 2 ,4  ,3 ) , DI V( 16, 5  J 

DATA    IWD1  ,R  WD/ 8, 9, 18, 18, 2  6,  2  8, 2, 2, 5,  5, 7,  8, 2.,  6.  ,10.  ,3.,  9.  ,15./ 

DATA  ALFA/1. DO, .9500/ 

DATA  AVA/-1.  ,-.50000000,0.0,.  250000000,1.    / 

AT  =  ALFAUNPI 

RATIO       =    .10000000 

IF(INP.NE.    2    .OR.  NUMSO   .LE.40)    GO    TO    10 

RATIO    =    FLOAT(NUMDE)/FLOAT{ NUMSO) 

CONTINUE 

DHMIN=HMIN/2. 

K*l 

CHECK  INPUT  PARAMETERS 

IFIMSTRT  .GT.  MAXO   .OR.   MSTRT  .LE.  0)  GO  TO  2000 

MO  =  MSTRT 

RNT=RINT(2)-RINT(1) 

IF(RNT  .LE.  O.DO)   GO  TO  2010 

IF(MAXO  .GT.  3   .OR.   MAXO  .LT.  1)   GO  TO  2300 

INITILIZE  PARTITION  VALUES  FOR  NEW  ORDER 


MAXN=MAXM( 

DO  30  1=1, 

NCK( I )=-l 

H=RNT/2. 

DP0S(M0,2, 

H=H/2. 

N=3 

SX(MO)=H 

NEQ(M0)=3 

NPT(l,l)=2 

NPT(1,2)=1 

NPT(1,3I=3 

DP0S(M0,2, 

DP0S(M0,2, 

DP0S(M0,1, 

DP0S(M0,2, 

DP0S(M0,1, 

DP0S(M0,3, 

DP0S(M0,1, 

DP0S(M0,2, 

DP0S(M0,3, 

DP0S(M0,3, 

IF( INP  .NE 


MO) 

MAXN 


1)=RINT(1)+H 


MAXN  + 
MAXN* 

2)  =  RI 
2)=RI 
1  )  =  DP 
2)=DP 
3)=DP 

3)  =  RI 
3)  =  RI 
1)  =  DP 
•  1) 


1)=RINT(1 I 
2)=RINT(2) 
NT(1  ) 
NTI1  )+H 
OS(MO,2,2) 
0S(M0,2,1) 
0S(M0,2,1) 
NT(2)-H 
NT  (2  ) 

OS(MO,2,3) 
GO  TO  1800 


STP139( 
STP140( 
STP14K 
,STP142( 
STP143C 
STP1440I 
STP1450( 
STP1460( 
STP1470(j 
STP14800 
STP1490C 
STP1500C 
STP1510C 
STP1520G 
STP15300 
STP1540C 
STP1550C 
STP1560C 
STP1570G 
STP1580C 
STP15900 
STP16000 
STP1610Q 
STP1620C 
STP1630C 
STP16400 
STP16500 
STP16600 
STP16700 
STP1680C 
STP16900 
STP17000 
STP17100 
STP17200 
STP17300 
STP17400 
STP17500 
STP1760Q 
STP17700 
STP178D0 
STP17900 
STP18000 
STP13100 
STP18200 
STP18300 
STP18400 
STP18500 
STP18600 
STP18700 
STP13800 
STP18900 
STP 19000 
STP19100 
STP19200 
STP19300 
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INITILIZES  ITEMS  FOR  THIS  PASS  OF  INTERVAL 


100 
110 


C 

c 
c 


120 


130 


210 


H=H/2. 

ILD»1 

INW»2 

JLD*1 

JNW=*1 

NMCK=0 

HH=H 

HHP=HH*. 941 59265 

HQ(1)*HH 

HQ(2I=HHP 

HSQ(1)=HH*HH 

HSQ(2»=HHP*HHP 

SETS  UP  THE  VALUES  FOR  ADDITIONAL  POINTS 

I=NPT(  ILD,JLD) 
IFINCKII I.GT.Ol  GO  TO   600 
FWD=DP0S(M0,2f I )+H 
BWD=DP0S(M0,2,I »-H 

CHECKS  TO  SEE  IF  THE  POINT  IS  ALREADY  IN  THE  PARTITION 

KA*-1 

DO  130  J=1,N 

IF<pABS<BWD-DP0S(M0,2,J    M.LT.DHMINJ    GO    TO      220 

IF(DABS(BWD-H-DP0S<M0,2,J    M.LT.DHMINI    KA  =  l 

CONTINUE 

IF(DABS(BWO-H-RINT(l»  KLT.DHMIN)    KA«1 

IF(KA    .LE.    01    GO    TO    2020 

CHECKS    TO    SEE    IF    MAXIMUM    NUMBER    OF    EQUATIONS    IS    EXCEEDED 

IF(N*1    #GT.    MAXNI       GO    TO    2D30 
N  =  N«-1 

lllMlllUi'rl  ?hEcke!tni  PARm,0N  P0,NI  ,NT0  PARTmoN 

DP0S(M0,2,N»=BWD 
NPT< INW,JNW|=N 
JNW=JNW*1 

DPOS(MOtl,N I ■DPOSIMOtl.il 

JC  =  N 

DP0S(M0,3fNI=DP0S(M0,2, I) 
JA=1 

IF(H.GT.HMAX)  GO  TO  560 

ERAR  =  EROR 

IF(  INP  .EO.  1»   GO  TO  300 

KLA=4 

XI=DP0S(M0,2,NI 
GO  TO  1850 
LUI  =  I 


STP19400 

STP19500 

STP19600 

STP19700 

STP19800 

STP19900 

STP20000 

STP20100 

STP20200 

STP20  30O 

STP20400 

STP20500 

STP20600 

STP20700 

STP20800 

STP20900 

STP21000 

STP21100 

STP21200 

STP21300 

STP21400 

STP21500 

STP21600 

STP21700 

STP21800 

STP21900 

STP22000 

STP22100 

STP22200 

STP22300 

STP2240G 

STP22500 

STP22600 

STP22700 

STP22800 

STP22900 

STP23000 

STP23100 

STP23200 

STP23300 

STP23400 

STP23500 

STP23600 

STP23700 

STP23800 

STP23900 

STP24000 

STP24100 

STP24200 

STP24300 

STP24400 

STP24500 

STP24600 

STP24700 

STP24800 
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L  (2  »=  NPTUNW,JNW-2I  STP24900 

L  (1 >  =  NPT( INW,JNW-3I  STP2500d 

1(51  =  NPT(ILD,JLD+1 )  STP25100 

IF(  JLD.GE.NEQ(MO)  I  L(5l=MAXN+2  STP25200 

IF1JLD.LE.2I  L(1I*MAXN+1  STP25300 

IF(JLD  .GT.  1)  GO  TO  400  STP2540Q 

L(2)  =  MAXN  +  1  STP2550(J 

L(1)  =  NPT<  ILD,JLDU)  STP25600J 

l(5l=NPT( ILD,JLD*2I  STP25700 

GO  TO  400  STP25800 

C  STP25900 

C      CHECKS  THE  PRESENT  PARTITION  POINT  AND  PUTS  POINT  TO  RIGHT  STP26000I 

C      INTO  THE  PARTITION  STP26100: 

C  STP26200 

220   NPT(INW,JNWI=I  STP26300 

JNW  =  JNW-H  STP26400 

REC=DP0S(M0,3,I»  STP26500 

DP0S(M0,lfI  )=BWD  STP26600 

0P0S{M0,3tI  )=FWO  STP26700 

IF(N+1.GT.MAXNI  GO  TO  2030  STP26800 

N=N*1  STP26900 

DP0S(M0,2fN»=FW0  STP27000 

NPTUNW,JNWl=N  STP27100 

JNW=JNW+1  STP27200 

DP0S(M0,l,NI=DP0S(M0,2tII  STP27  300 

DP0S(M0,3tN)=REC  STP27400 

JC=N  STP27500 

JA  =  3  STP27600 

ERAR  =  EROR  STP27700 

IFUNP  .EQ.  II   GO  TO  300  STP27800 

KLA=5  STP27900 

XI=DP0S(M0,2,NI  STP28000 

GO  TO  1850  STP28100 

230   L(1)=NPT( INW,JNW-3I  STP28200 

L(2)=NPT( INW.JNW-21  STP28300 

L(4I=NPT<  ILD,JLD+1>  STP284CG1 

L(5)  =  NPTULDt  JLD+2)  STP28500 

IF(  JLD.GE.NEQ(M0)-1  I  l(5l=MAXN*2  STP28600' 

IF( JLD.LT.NEQ<M0I I  GO  TO  430  STP28700; 

L(4)  =  MAXN*2  STP28800 

L(5)=NPT< INW,JNW-4I  STP28900 

GO  TO  400  STP29000 

C  STP29130 

C      CHECKS  THE  POINT  TO  THE  RIGHT  OF  THE  OLD  PARTITION  POINT.  STP29200 

C  STP29300 

240   JA=2  STP29400 

JC=I  STP29500 

ERAR  =  EROR  STP2960C 

IFtINP.EQ.il  GO  TO  300  STP29700 

XI=DP0S(M0,2»JC»  STP29800 

L(5>  =  NPT(  ILD,JLD-HI  STP29900 

L (2>=NPT( INW.JNW-31  STP3D000 

L(4»=NPT( INWfJNW-ll  STP30100 

L(l  »  =  NPT(  INW,JNW-4)  STP30200 

IFULD.EQ.ll  L(1I  =  MAXN+1  STP30300 
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IF(JLD.GE.NEQ<MO)l  L(5I»MAXN*2  „„r,nn 

c      G0  T°   «°°  liplltlo 

C     FINOS  FUNCTIONAL  VALUES  FOR  INITIAL  PARTITION  STP30700 

C300   DO  310  JB-1,3  |JJ^JJ 

VI  =  DP0S(M0,2,JU*AVA(  JB>*HH  STP31000 

IF(DABS(V1-RINT(1H.LT.  1.1*HH  .OR.  DABS  ( Vl-RINT  (2)1    .  LT.       STP3110Q 

1   1.1*HHJ    ERAR  =  .5  *  ER3R  zll.Z.iz.Z. 


STP31200 


1.1*HH)    ERAR  ■  .5  *  EROR 

CALL  INIT(A(1,1,JB),M0,V1J  STP^lino 

310   CALL  INIT(A(1,2,JB>,M0,DP0S(M0,2,JC>*AVA(JB»*HHP»  STP31400 

GO  TO  (320, 350, 380), MO  STP31500 

C      EVALUATES  USING  LAGRANGIAN  -  INITIAL  STP31700 

C320   JB'°  STP319001 

330   JB  .  JB  ♦  I  iTP320oS 

IF(IDER(1,1I.EQ.0)  GO  TO  340  STP3710O 

AA(ll  =  (A(l,JB,5l-A(l,JB,in/HQ(JB»  STP 32200 

,JCr!f,/!inl,,o?,^^^?<?J0Q0<tA(l,JB,3,  +  16'*A(1*JB»2,-6'000000D*A«l»JSTP32300 

1  a  »  1  I  I  /  {  n\J  i  J  B  i  *H  5Q  I  J  B  II  ^TPUtnn 

AE(1)=AA<H*AC/6.0000000  *HSQ(JB>  STP^nn 

340   IF(IDER(2,1>.EQ.0>  GO  TO  510  STP32600 

AA(2I  =  (A<1,JB,5|-2.*AU,JB,3)*A(1,JB,1U/HSQ(  JB  »  STP 32 700 

AC=( 10. 666666666666667' All, J  8,5  I -136. 53333333 3333* A (1,JB, 4 » +192. 00 STP32800 

1000*A(1,JB,3J-85.3333333333333*A<1,JB,2)*19.200000*A{1,JB,1I)/<HSQSTP32900 

2IJBI*HSQ(  JBH  STP3300O 

AE(2I=AA(2»*AC/12. OOOOOO  *HSQUBJ  STP33100 

c      G°  T°  5l°  STP33200 

C      EVALUATES  USING  HERMITE     -  INITIAL  STP33400 

C  isn    in  -  n  STP33500 

360   JB  -  m  ♦  l  STP33600 

Al'll-ltl    \«    *.  STP33700 

AA( 1)-A(2,JB,3»  cTP,,ftM 

AE(1I=AA(1»  ilKjjBOO 

IF<IDER(2,2I.EQ.0I  GO  TO  370  STP34000 

1AA(,f,;,?*!!A,1,JR,5,"2,*A(1,JB,3,*A(1'JB'1M/HSQ(JB'-(A<2.JB,5)-A(STP34100 
J.l,JB,ll!  /(2.*HQ(JB))  STP3420Q 

,A^.!T?Mn3o33^3333333*A{1,JB,5'"8640'0000*A(1»JB'3»*8533.33333333STP34300 

13  33*A(l,JB,2»*72G.OOOOO*A<l,JB,in/HSQ(JBI*(160.00COG*A<2,JB,5>*28STP34400 

280.0JCC*A(2,JB,3l*2560.00O0*A(2,JB,2N  /HQ(JB»  STP34500 

AE<2>=AA(2)*AC/32  0.  00000  stpiZa™ 

37C   IF(IDER(3,2I.E0.CI  GO  TO  510  STP34700 

/A<3j=<(A(ltJB.5J-A(l,JB,in*7.5000000/HQ<JB»-1.5000000*(A(2,JB,5ISTP34800 
l*8.*A(2,JB,3»*A(2,JB,im/HSQUB»  STP34900 

AC= (( 352 80. 000* A( 1 , JB, 1 1*47786. 6666666667* A  II, J  8,21-80640. 000*A( 1 , STP35000 
J;)8!  3>  "242  6.  66666666667*  A(1,JB,5>  I /HQ<  JB » * ( 5040. O0OC*A (2 , JB , 1 1 +3584STP35100 
10.000*A(2,JB,2I*20160.000*A(2,JB,3I*560.00000*A(2,JB,5M>/HSO(JB»  STP35200 

rnVn    ^a**3'  *  AC/84O.OOO000  STP35300 

c      G°  T°  51°  STP35400 

C      EVALUATES  USING  TRIMITE     -  INITIAL  STP35600 

380   JB  -  O  STP35700 

380   JB  "  °  STP35800 


C 
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390   JB  »  JB  ♦  1  STP3590 

AA(1)=A(2,JB,3I  STP3600 

AE(1I=AA(1I  STP3610 

AA(2I=A(3,JB,3I  STP3620 

AE(2I*AA(2I  STP3630 

IF(  IDER(3,3).EQ.O)  GO  TO  395  STP3640 

AA(3»»ll3.1250L0/HSQ(JB)*CA(l,JB,5»-A(l.JB,l) 1-1. /HQ( JB I  *  (4.  125000STP3650 

10*A(2,JB,5»+18.00  0000*A(2,JB,3H-4.1250000*A(2,JB,1)  K.  375  00000*  ( A  (  STP3660 

23t  JBt  5I-A(3f  JBtim/HQ(  JB)  STP3670 

AC=(  (6441 12  00.  *A(  1,  JB,H- 263782400.  *A(  I,  JB,  2) +199584000.  *A(1,JB,  3)  STP3680 

l-212  8uC.DO*A(l,JB,5  I  » /HSQ( J B » ♦( 1065960C. * A( 2, JB, 1 1-344C6400. *A( 2, JSTP3690 

2B, 2 ) -43545600. *A( 2, JB, 3  J  +  25200. 000* A( 2 , JB ,5  I » /HQ( JB )♦ ( 498960. 00* A ( STP3700 

33, J B,l) -9461 760, 0*A  (  3,  JB,  2) +3265920.  0*A(  3  ,  JB,  3H-1680.  0000*A  (3  ,  JB,  5STP3710 

4)>>/HQ(JB)  STP3720 

AE( 3»=AA(3»+AC/60480.000  STP3730 

395   IF( IDER(4,3».EQ.0»  GO  TO  510  STP3740 

AA(4l=(72.000000/HSQ(JB)*(A(l,JB,5>-2.*AU,JB,3>+A(l,JB,l»)-19.500STP3750 

10000/HQ(JBI  STP3760 

1*(A(2,JB,5)-A(2,JB,1)»+1. 5000000* (A(3,JB, 51-24. OOOOOQ*A( 3 ,JB, 3 J+A( STP3770 

33,JB,1  m/HSQ(JB>  STP3780 

AC=(( 1146880000. *A( 1, JB , 2 »-l 5240960C. *A( 1 , JB, 1 » -1001548800, *A ( 1 , J BSTP3790 

1,3»+7078400.0*A(1,J8,5I J/HSQl JB >+( 206438400.* A( 2, J B, 2 1-22680000. *ASTP380Q 

3  (  2, JB, 11+2  32243200. *A(2 , J B, 3l-204963D.0*A< 2, JB,5 I » /HQ( JB )+ ( 344C64CSTP3810 

30.*A(3,JB,2»-907200.00*A(3,JB,1>-21772800.*A(3,JB,3»+168000.00*A(3STP3820 

4,JB,5)n/HSU(JB)  STP3830 

AE(M=AA(4)+AC/151200.00  STP3840 

GO  TO  510  STP3850 

C  STP3860 

C      PREPARE  FOR  EVALUATION  OF  THE  DERIVATIVE  STP3870 

C  STP3880I 

400   L(3»=JC  STP3890 

DO  410  IB=1,5  STP3900i 

LA=L(IB)  STP3910 

XT(  IB)=DP0S(M0,2,LAI  STP3920 

DO  41C  IA=1,M0  STP3930 

410   YT(  IA,IB»=YTEMP(IA,LA»  STP3940 

CALL    DERIV(XT,YT, IDERd ,M0) , MO, 0, 0 ,D       , MOLD,X I , GAM, HQ( 1 >>  STP3950 

GO    TO    (420,440,4601 , MO  STP3960 

C  STP3970' 

C  EVALUATES    USING    LAGRANGIAN   -    MODIFYING  STP3980 

C  STP3990 

420   IF( IDERd, D.EQ.OI  GO  TO  430  STP4000 

A A(  1  I  =  (YT(1,4)-YT( 1,2)1/(2.  0*HQ(1)>  STP4010' 

AE(  1)=AA(1)+DVY<1  I  STP402CI 

430   IF( IDER(2,l).EQ.O)  GO  TO  500  STP4030I 

AA(2»=(YT(1  ,4)-2.*YT(l,  3)+YT(l  ,2))/HSQ(l)  STP4040( 

AE(2J=AA(2)+DVY(2  »  STP4050I 

GO  TO  500  STP4060 

C  STP4070I 

C      EVALUATES  USING  HERMITE     -  MODIFYING  STP4080I 

C  STP4D90( 

440   AA(1)=YT(2,3»  STP4100I 

AE(1)=AA(1>  STP4U0 

IF(  IDER(2,2I.EQ.0I  GO  TO  450  STP4120* 

AA(2)  =  2.*(YT(1  ,4I-2.*YT(1  ,3)+YT(l  ,2 )) /HSQ ( 1 »-( YT(  2  ,4I-YT(2  ,  STP4130I 
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450 


460 


470 


500 
510 

520 


530 
540 


1211/(2. 0*HQ(1N 

AE(2»=AA(2»+DVY(2) 

IF( IDER(3,2).EQ.O)    GO   TO    500 

l"i?(,2,:3IUT;2,;II!}/Hsi:n•50oooooo/Ha(l,-l■5o3oool>•|YT,2 

AE(31=AA(3»*DVY(3I 
GO    TO    500 

EVALUATES    USING   TRIMITE  -    MODIFYING 

AA(1)=YT(2,3I 

AE(1I*AA(1) 

AA(2)=YT(3,3I 

AE(2)=AA(2> 

IF<  IDER(3,3).EQ.OJ    GO    TO    470 

AA(3»  =(  13.1250uC/HS0(l»«-{  YT(1 
2YT(2,2n+18.CC0000*YT(2,3M/ 
1HQ(1»+. 37500000*1 YT(3    ,4>-YT(3 

AE(3>=AA(3I*DVY(3» 

IF(IDER(4,3».EQ.0I    GO    TO    500 

AA(4>=(72.0000GO/HSQ(l»*(YT(l  ,4)-2.*YT(l 
1HQ(1>»(YT(2  ,4)-YT(2  ,2 ) I *1. 5000000* ( YT( 3 
2TC3    ,2MI/HSQ(ll 

AE(4)=AA(4»+DVY(4I 

MAKES    USE    OF    THE    USER'S    EQUATION 

CALL    UFUN(M0,DP0S(M0,2,JC),YT(1,3»,T,AA,AD,1I 

CALL    UFUN(M0,0P0SCM0,2,JCI,YTC1,3J,T,AE,AF,1I 
GO    TO    520 

CALL    UFUN(M0,DP0S(M0,2,JC»,A(1, JB,3>,0. ,AA,AD,1I 
CALL    UFUN(MO,DPOS(MO,2,JC),A(l,JB,3»,0.,AE,AF,l| 

CHECKS    THE    ERROR 

DO    540    IL    =*    1,M0 

AF( IL)    =    AF{ IL I    -    ADCILI 

IF(0ABS(AD(  ID  »    .GT.    1.D0I 

AD(  IL)    =    AD(ID/AD(  ID 

IF(DABS(AF(IL»)     .GT.    ERARI 

CONTINUE 

IF{  INP    .EQ.    2      .OR.       JB    .GE.    2) 

GO    TO    (330,360,390), MO 


,4>-YT(l    ,2))-(4.1250000*(YT(2 
,2>))/HQ(l> 


,3»+YT(l    ,211-19.5 
,4)-24.  000000*YT(3 


GO   TO    533 

GO   TO    560 


GO    TO    550 


!??iniTES    ™AT    ™E    ERR0R    IS    ALL    RIGHT,    SO    THAT    THE    POINT    IS 
IGNORED    ON    FUTURE    PASSES 


550 
560 


NCK(JC)    »    1 

NMCK    =    NMCK    ♦    1 

GO    TO    (220,610, 2401, JA 

CHECKS    FOR    COMPLETION    OF    PASSES. 
FOR    NEW    PASS 


600   NPT(INW,JNWI»I 


IF  NOT  DONE,  PREPARES 


STP41400 
STP41500 
STP41600 
,41+8  STP41700 
STP41800 
STP41900 
STP42000 
STP42100 
STP42200 
STP42300 
STP42400 
STP42500 
STP42600 
STP42700 
STP42800 
»4)  +   STP42900 
STP43000 
STP43100 
STP43200 
STP43300 
000GD/STP4340O 
,3I*YSTP43500 
STP43600 
STP43700 
STP43800 
STP43900 
STP44000 
STP44100 
STP44200 
STP44300 
STP44400 
STP44500 
STP44600 
STP44700 
STP44800 
STP44900 
STP45000 
STP45100 
STP45200 
STP45300 
STP45400 
STP45500 
STP456T)0 
STP45700 
STP45800 
STP45900 
STP46000 
STP46100 
STP46200 
STP46  300 
STP46400 
STP46500 
STP46600 
STP46700 
STP46800 
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JNW*JNW*1 

NMCK=NMCK*1 
610       IF( JLD.GE.NEQ(MO) I    GOTO      620 

JLD=JLD+1 

GO    TO    123 
620      IF(NMCK    .GE.    Nl    GO    TO   790 

ILD=3-IL0 

INW=3-INW 

NEQ(MO>»N 

H=H/2. 

IF(H    .LT.    HMIN)       GO    TO   2090 

SXl  MO)»H 

GO    TO    110 

: 

C  PREPARE    THE    WORK    ESTIMATE 

: 

700   IF(NEQ(MOI  .GT.  0  1   GO  TO  710 

WVARJMOI  =  1.0050 

GO  TO  720 
710   WVTP  =  IDERIM0,M0)*IWD1(M0,U  ♦  IDER( M0+1,M0I *IWD1 (M0,2) 

TMO  =  2*M0 

WVTP    =    WVTP/TMO 

WVTR    =    IDERJM0,M0I*IWD1(M0,3>    ♦    I  OER  (  MO  +  1  ,M0 )  *I  WOK  MO  f4 1 

RN    =    N 

WVAR{MO>  =  RN*(WVTP  ♦  RAT  10* i WVTR*RN+RWD( M0,1 » )  ♦  RWD<M0,2») 
720   IF(MO  .EQ.  MAXO )  GO  TO  730 

MO  =  MO  ♦  1 

k=i 

GO  TO   20 
:      DECIDE  WHICH  ORDER  TO  USE 

730   MO  =  MSTRT 

IF(MSTRT    .EQ.    MAXOJ    GO   TO   740 

IF(WVAR(M0)*AT    .GT.    WVAR(M0  +  1M     MO   »    MO    ♦    1 

IF( MAXO-MSTRT    . EQ.     1)    GO    TO    740 

IF(MO    .EQ.    DAT    =    AT   *AT 

IF( WVAR(MO»*AT  .GT.  WVAR(3>I  MO  =  3 
74-0   MONW  =  MO 

IF(WVAR(MO»  .GE.  1.0D45)   GO  TO  2040 

RETURN 

CHECKS  THE  PARTITION  TO  INSURE  SMALLEST  AVAILABLE  STEP 
SIZE  IS  USED 

790   JNW  =  1 

NEQ(MO)=N 
800   JA=NPT( INW,JNW» 

IF<  JNW.EQ.l)  GO  TO  805 

IF( JNW.NE.NI  GO  TO  810 
H=RINT(2»-DPOS(MO,2,JAI 

JLD=-1 

GO  TO  830 

805   H=DP0S(M0,2,JAI-RINT(1I 
JLD=1 


STP4690I 
STP4700< 
STP4710I 
STP4720I 
STP4730< 
STP4740I 
STP4750< 
STP4760( 
STP4770C 
STP4780< 
STP4790( 
STP4800C 
STP4810C 
STP4820C 
STP4830C 
STP4840C 
STP4850C 
STP4860C 
STP4870C 
STP4880C 
STP4890C 
STP4900C 
STP4910C 
STP4920C 
STP4930C 
STP4940Q 
STP49500 
STP4960C 
STP49700 
STP4980Q 
STP49900 
STP50000 
STP50100 
STP5D200 
STP5O3O0 
STP50400 
STP50500 
STP5C60C' 
STP5070C 
STP538D0 
STP50900 
STP51000 
STP51100 
STP51200 
STP51300 
STP51400 
STP51500 
STP51600 
STP51700 
STP51800 
STP51900 
STP52000 
STP52100 
STP52200 
STP52300 
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810 


820 
830 

840 
850 

860 


GO    TO    830 
JB»NPT(INW,JNW-1» 
JC=NPT< INW,JNW*1) 
JLD  =  1 

VA=DPOS(MO,2,JA)-DPOS<MO,2, JB) 
IF{  VA.LE.O.  I    GO    TO  2050 

AB»DP0S(M0,2, JCI-0P0S(M0,2t JAI 

IF(  AB.LE.O. I    GO    TO    2350 

IF(DABS(VA-AB).LE.DHMIN»    GO    TO    850 

IF(VA.GE.AB|    GO    TO      820 

JLD=-1 

H=*AB 

GO    TO    830 

JLD=i 

H=VA 

DO    840    ILD»1,N 

JB= JNW*JLD*ILD 

IF(  JB.LE.O.OR.JB.GT.NI    GO    TO    2060 
JB=NPT( INW, JBI 


JoN??NUEDABS(DPOS<MO,2,JB'"DPOStM°'2'JAn-H,-LT-HMIN/2-»    GO    TO 


DFLOAT(  JLD*(ILDH 


C 

GO    TO    2060 

DP0S(M0,1,JA) 

ILD»1 

JNW=JNW+1 

IF(JNW.LE.N)    GO    TO      800 

DO    86C    1  =  1, N 

NPTER(MO,II*NPT( INW.IJ 

GO    TO    700 

CHECKS  THE  OLD  PARTITION  FOR  VALUES. 
IT  INTERPOLATES 


IF  NOT  PRESENT, 


180C 


IF(MO  ,GT.  MOLD!  GO  TO  2070 

JC=MAXN+1 

N3»l 

KLA  =  7 

GO  TO  1930 
1810  JC=MAXN*2 

N3=NUM+2 

KLA  *  6 

GO  TO  1930 
1820  KLA*2 

XI*. 25 

JC=»2 

GO  TO  1850 
1830  KLA=3 

XI=. 7500000000 

JC  =  3 

GO  TO  1850 


1840 
18  50 


KLA  =  1 

JC=1 

XI».5 

Nl  =  l 

N2*NUM 


♦  2 


STP52400 
STP52500 
STP52600 
STP52700 
STP52800 
STP52900 
STP53000 
STP53100 
STP53200 
STP53300 
STP53400 
STP53500 
STP53600 
STP53700 
STP53800 
STP53900 
STP54000 
STP54100 
STP54200 
850STP54300 
STP54400 
STP54500 
STP54600 
STP54700 
STP54800 
STP54900 
STP55000 
STP55100 
STP55200 
STP55300 
STP55400 
STP55500 
STP55600 
STP55700 
STP55800 
STP55900 
STP56000 
STP56100 
STP56200 
STP56300 
STP56400 
STP56500 
STP56600 
STP56700 
STP56800 
STP56900 
STP57000 
STP57100 
STP57200 
STP57300 
STP57400 
STP57500 
STP57600 
STP57700 
STP57800 


.9.2. 


1860    N3»(Nl*N2)/2 

IF(  N3.EQ.N1  » 
IF  (  X(N3I  - 
N1«N3 


GO   TO    1890 

XI)    1870*1930*1880 


1870 
1880 
1890 


1900 
1910 


1920 

1930 
1940 
195C 
2000 

2010 

2020 

C 
C 
C 

2030 

2040 
2050 
2060 
2C70 
2080 
2090 


GO    TO    1860 
N2-N3 

GO  TO  1860 
IF(N3  «LT. 
ITST  *  NUM 
IF(N3  *EQ. 
IF(N3  .EQ. 
00  1910  111 
112    »    ITST 


1    «0R.    N3    .GT.    NUMU)    GO    TO    2080 

-2 

1    I    ITST   *  0 

NUM    ♦    11     ITST    «    NUM   -    2 

»    1,4 
♦    111 


DO    1900    113=1, MOLD 

YT(  I13,I11»*Y(1,I13,I12> 

XT( I11)=X(I12I 

CALL    DERIV(XT,YT,IDER(1,M0),M0,0,1,       0. MOLD, XI, GAM, HQ(1 II 

DO    1920    113    =    1,M0 

YTEMP(I13,JC)»GAM( 1 13  I 

GO    TO    1950 

DO    1940    113    *    lfMO 

YTEMP(U3,JC»*Y(1,I13,N3) 

GO  TO  (1820, 1830,100, 210,230, 1840, 18101 *KLA 

K  »  -201 

RETURN 

K  »  -202 

RETURN 

K  =  -205 

RETURN 

RETURN  COOES  ARE  GENERATED 

K  =  -206 

NEQ(MOI»-3 

GO  TO  700 

K  *  -209 

RETURN 

K   =  -208 

RETURN 

K  =  -207 

RETURN 

K  =  -203 

RETURN 

K  =  -204 

RETURN 

NEQ(MO)  =  -213 

GO  TO  7C0 

END 


STP57900 
STP58000 
STP58100 
STP58200 
STP58300 
STP58400 
STP58500 
STP58600 
STP58700 
STP58800 
STP58900 
STP59000 
STP59100 
STP59200 
STP59300 
STP59400 
STP59500 
STP59600 
STP59700 
STP59800 
STP59900 
STP60000 
STP60100 
STP60200 
STP6D300 
STP60400 
STP60  500 
STP60600 
STP60700 
STP60800 
STP60900 
STP61000 
STP61100 
STP61200 
STP61300 
STP61400 
STP61500 
STP61600 
STP61700 
STP61800 
STP61900 
STP62000 
STP62100 
STP62200 
STP62300 
STP62400 
STP62  500 
STP62600 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


SUBROUTINE    DERIVU,  FN,  IDER,MO,KB, KF,KG,  MOLD,  XI. GAM, HI 

THIS    PROGRAM    IS    USED    TO   FORM    INTERPOLATION    POLYNnMTAK    cno 
LAGRANGIAN,    HERMITE,    AND    TRIMITE    INTERPOUtJS    POLWOHuS 
FOR    FUNCTIONAL    VALUES    AS    WILL    AS    SEVERAL    DERIVATIVES         li 
IS    DIVIDED    INTO    TWO    SECTIONS.       ONE    SECTION    CALCULATES    ill 
INTERPOLATING    POLYNOMIALS    NEEDED    IN   EST    MATING    ERROR I         THf 
OTHER    SECTION    IS    USED    TO    CALCULATE    NEW    VALUES    FOR    fSnCTiJnI 
DESlRll0ClAJED    DERIVATIVES    WHEN   NEW    PARTITION    POINTS    ARE 


THE    VARIABLES 

X 

Y 

IDER 

MO 

MOLD 

ALPHA 

DIV 


COEF 
KB 


IN    THIS    SUBPROGRAM    ARE 
CONTAINS    INTERPOLATING    POINTS 
CONTAINS    FUNCTION    AND    DERIVATIVE    VALUES 
TELLS    WHICH    DERIVATIVES    ARE    REQUIRED 
ORDER    OF    METHOD    (POSSIBLY    NEW) 
OLD    MO 
CONTAINS 
FOR    EACH 
CONTAINS 
(1,1) 
(2,11 
(3,1) 
C4,IJ 
(5,1) 
(6,  I) 
(7,1) 
(8,11 
(9,  II 
(10,1  I 
(11,1  I 
(12,1) 
(13,1) 
(14,1  ) 
(15,1) 
(16,1  ) 
CONTAINS 


THE    FACTORS      X-A(I)    AND    A(J)-A(II 

POINT    AND    COEFFICIENT 

VARIOUS    FUNCTIONAL    VALUES    FOR    EACH    POINT 
L"MXI 
LMX) 

LMIMX(II) 
L'MXI 

L'MIIIXIIII 
DENOMINATOR    OF    L(X) 
A* 

(A    BAR)* 
(A    DOUBLE    BAR)* 

A*  • 

(A    BAR)* • 

DOUBLE    BAR)** 


(A 
A 
A 
A 


BAR 

DOUBLE    BAR 
L(X) 
THE    COEFFICIENTS 


THOSE    FROM 
7   TO      15 
ARE    ALSO 
USED    TO 
STORE 
TEMPORARY 
RESULTS    IN 
INTERPOLATION 


OF 


THE  RESULTING 
OTHERWISE 


KG 


KF 

XI 

GAM 

DVY 

H 


1     IF    COEFFICIENTS    DESIRED,    0 
INTERPOLATING    POLYNOMIALS 
1-INTERPOLATE       0-ERROR    ESTIMATE 
1-AT    PARTITION    POINT         0-NOT    AT 
EVALUATION    POINT 
RESULTS    OF    INTERPOLATION 
RESULTS    OF    ERROR    ESTIMATES 
DISTANCE    USED    IN    ANALOF.     HSQ,HSQQ ,&HSQC 
rur    occr    nc    -<H**2),H***,H»*6    RESPECTIVELY 
THE    REST    OF    THE    VARIABLES    ARE    POINTERS    AND   COUNTERS 

IMPLICIT  REAL  *  8  (A-H.Q-Z) 
DIMENSION  X(1|,FN(3,5>,IDER(4|,GAM(3I 

HSQ  =  -  H  ♦  H 


ARE 


DRV10000 

DRV10100 

DRV10200 

DRV10300 

DRV104OO 

DRV10500 

DRV10600 

DRV10700 

DRV10800 

DRV1C930 

DRV11000 

DRV11100 

DRV11200 

DRV11300 

DRV11400 

DRV11500 

DRV11600 

DRV11700 

DRV1180O 

DRV11900 

DRV12000 

DRV12100 

DRV12200 

DRV12300 

DRV12400 

DRV12500 

DRV12600 

DRV12700 

DRV12800 

DRV12900 

DRV13000 

DRV13100 

DRV1320O 

DRV1330O 

DRV13400 

DRV13500 

DRV13600 

DRV13700 

DRV13800 

DRV13900 

DRV14000 

DRV14100 

DRV14200 

DRV14300 

DRV14400 

DRV14500 

DRV14600 

DRV14700 

DRV14800 

DRVK900 

DRV15000 

DRV15100 

DRV15200 


9k 


10 

c 
c 
c 


c 
c 
c 

20 


C 

c 
c 

30 


40 


C 
C 
C 


45 


50 

C 
C 
C 


IF(M0  .LT.  2  »  GO  TO  5 

HSOQ*HSQ*HSQ 

IFCMO  .LT,  3  I  GO  TO  5 

HSQC»HSQQ*HSQ 

TEST  IF  DERIVATIVES  ARE  USED 

KA  *  0 

DO  10  IA»1,4 

DVY( IA»=0. 

KA=KA*IDER{  IA  I 

IF(KA.  EQ.O  )  RETURN 

SETS  FOR  HIGHEST  LAGRANGIAN  DERIVATIVE 

IF(M0.NE.1.0R.IDER<2).NE.l)  GO  TO  20 

Nl»5 

K=2 

VA*12. 

V=24. 

GO  TO  30 

SETS  FOR  THE  LOWER  LAGRANGIAN  DERIVATIVE 

IF(M0.EQ.1.AND.IDERU).EQ.  0)  RETURN 

K»l 

Nl=4 

VA=6. 

V*6. 

OBTAINS    DENOMINATOR    OF   LAGRANGIAN    POLYNOMIALS 

DO    40    IA=1,N1 

DIV(6,IA»-1. 

DO    40    IB=1,N1 

IF{  IA.EQ.  181    GO    TO   40 

DIV(6,IA»=DIV(6,IAI*(X(IAI-X(IBI» 

CONTINUE 

IF(KF.NE.l)    GO    TO    45 

KSAVE=KA 

GO    TO    70 

THE    FOLLOWING   FORMS    L^'MXI 

DO    50    IA=1,N1 
DIV(1,IAJ=V/DIV(6,IA) 
IF(M0,NE.1I    GO    TO   70 

CALCULATES  ERROR  ESTIMATES  FOR  LAGRANGIAN 

DVY(K>*0. 

KA  =  4*(K-1) 

DO  60  IA*1,N1 

IF(KB  .EO.  0  I  GO  TO  58 

IB  =  KA  ♦  IA 


DRV1530I 

DRV1540I 

DRV1550I 

DRV1560! 

DRV1570* 

DRV1580< 

DRV1590<; 

DRV1600I 

DRV1610< 

DRV16204 

DRV1630I 

DRV1640<! 

DRV1650( 

DRV1660C 

DRV1670<\ 

DRV1680(] 

DRV1690C! 

DRV1700C 

DRV1710C 

DRV1720CJ 

DRV1730C 

DRV1740C 

DRV1750( 

DRV1760C 

DRV1770C 

DRV1780C 

DRV1790C 

DRV1800C 

DRV1810C 

DRV1820C 

DRV1830C 

DRV1840C 

DRV1850C 

DRV1860C 

DRV1870C 

DRV1880C 

DRV1890C 

DRV1900C 

DRV1910C 

DRV1920C 

DRV19300 

DRV194O0 

DRV19500 

DRV19600 

DRV19700 

DRV19800 

DRV19900 

DRV20000 

DRV20100 

DRV20200 

DRV20300 

DRV20400 

DRV20500 

DRV20600 

DRV20700 
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C 

c 
c 

70 


C0EF(  IB,1,1I«DIVI1,IAI/VA      *HSQ 
GO    TO    60 
58      DVYIKI    -    DVYfKI    ♦    DIV( 1 , I AI*FN( 1, IA) 
60  CONTINUE 

IF(KB    .EQ.l I    GO   TO   62 
DVY(K)»DVY(K»/VA  +HSQ 

62      IF(K    .EQ.    1)    RETURN 
GO    TO    20 

CALCULATES    THE    FACTORS    USED    IN   THE    NUMERATOR 

DO    80    IA*l,Nl 

KA=l 

DO    80    IB=1,N1 

IFI  IA.EQ.  IB)    GO    TO    80 

ALPHAI1,  IA,KAI=XI-XUBI 

ALPHA12,IA,KA1*X<  IAI-XIIBI 

KA=KA  +  i 

CONTINUE 

IFIKF.EQ.ll     KA=KSAVE 

IFIKF.EQ.ll    GO    TO    I  270, 90,120 1 , MOLD 

KA*2 

CALCULATES      LMXI      OR      LMIMXIIII 

DO    100    IA=l,Nl 

DI  V  (KA*1,IA  I  =  (  ALPHA  (KA,IA,U»(ALPHA(KA,IA,  21*  ALPHA! 
1IKA,IA,2I*ALPHAIKA,  I A ,3 ) ) /DI V (6 , I  A  I 
IFIKF.EQ. II    GO    TO    280 
IFIKA.EQ.ll    GO    TO    220 
IFIM0.LE.2I    GO    TO    140 

CALCULATES      L'MXI       OR      L'MIMXIIli 

KA*1 

DO    130    IA=l,Nl 

D IV (<A*3,IA)  =  (2,*( ALPHA (KA, I A.l !♦ ALPHAI KA, I  A, 21 +ALP 
1IVI6,IA) 
IFIKF.EQ.ll    GO   TO    90 
IFIKA. EQ. 21    GO    TO    200 
IF( MO. EQ. 2  I    GO    TO    180 
KA»2 
GO    TO    120 

CALCULATES    A«     FOR    HERHITE 

DO  150  IA»l,Nl 

DIVI7,IAI=-2.*DIVI3,IA| 
DIVI8,IA)*1. 

IFI  IDERI2I.NE.0I  GO  TO  110 

CALCULATES  THIRD  DERIVATIVE  FOR  HERMITE 

IFI  IDERI3I.EQ.0I  RETURN 
DO  170  IA*1,N1 


80 


C 

C 

C 

90 

103 


C 

c 
c 
no 

120 
130 


DRV20800 
DRV20900 
DRV21000 
DRV21100 
DRV21200 
DRV21300 
DRV21400 
DRV21500 
DRV21600 
DRV21700 
DRV21800 
DRV21900 
DRV22000 
DRV22100 
DRV22200 
DRV22300 
DRV22400 
DRV22500 
DRV22600 
DRV22700 
DRV22800 
DRV22900 
DRV23000 
0RV23100 
DRV23200 
DRV23300 
KA,IA,3I I+ALPHADRV23430 
DRV23500 
DRV23600 
DRV23700 
DRV23800 
DRV23900 
DRV24000 
DRV24130 
DRV24200 
DRV24300 
HAIKA,IA,3II)/D0RV24400 
DRV24500 
DRV24600 
DRV24700 
DRV24800 
DRV24900 
DRV25000 
DRV25100 
DRV25200 
DRV25300 
DRV25400 
DRV25500 
DRV25600 
DRV25700 
DRV258C0 
DRV25900 
DRV26000 
DRV26100 
DRV26200 
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171 
170 


C 
C 
C 
180 


C 
C 
C 
200 


210 


C 
C 

c 

220 


V-140. 00000+0  I V(l , IA»*DIV(1,IA| 

COEF(  IA,2,1I«DIV(7,IAI*V 

COEF(IA,2,2l=DIV<8, IAI*V 

IF(KB.EQ.O)    GO   TO    171 

COEFdA, 2, II    *    COEFdA, 2,  11/840.000000   *    HSQQ 

COEFdA, 2, 2)    *   COEFdA,  2, 21 /84C. COOOOO    *    HSQQ 

GO    TO    170 

DVYm  =  DVY(3l*FNd,IAI*COEFdA,2,ll  +  FN{2,IAI*COEFdA,2,2) 

CONTINUE 

IFIKB    .EQ.    II    RETURN 

DVY(3I    *    DVY(3I    /    840.000030    *    HSQQ 

RETURN 

CALCULATE  A   FOR  HERMITE 


*DIV(3,IA| 


DO    190    IA*1,N1 
DIVd3,IAI=l-2.*(XI-XdA) I 
DIVd4,IAI=<XI-X< I A } I 
V*20.000000*DIVd,IAI*DIVd»IAI 
VA»120.03000*DIV(4,IAI*DIV(1,IAI 

;  CALCULATES    SECOND    DERIVATIVE    FOR    HERMITE 

COEF(  IA,1,1I»V*DIV(  13,IAI+VA*DIV(7,IAI 

COEFdA,  1,2  I  =V*DIV( 14, 1  A I +VA*DI V( 8, I  A) 

IF(KB.EQ.O»    GO    TO    191 

COEFJ I  A, 1,1  I =COEF( I  A, 1,11/360. 00000  *HSQQ 

COEFdA,l,2)=COEFdA,l,2l/360.00000     *HSQQ 

GO    TO    190 
191      DVY(2l=DVY(2l+FNd, IA»*COEF(IA, 1 ,  ll  +  FNl  2,  IA  l*COEFU  A,  1,21 
190      CONTINUE 

IF    (    KB    .EQ.lt    GO   TO    160 

DVY(2)=DVY(2l/360.00000  *HSQQ 

GO   TO    160 


CALCULATES    A«,A»»    FOR   TRIMITE 

DO    210    IA=1,N1 

DIV(ll,IA)=-6.C00000n*0IV(3,IAI 

DIV(10,IA)=-2.*DIV( 11,IAI*DIV(3,IAI-3.00000000*DIV(5,  IAI 

DIVd2,IA)=l. 

DIV(7,IA»=DIVdO,  IAI*(XI-X    dAI)  -3. 0000000*01 V( 3, IA I 

DIV(8,IA»=1.  +  (XI-X(  IAI  I     tDIVtllf  IA! 

DIV(9,IA)=XI-X( IAI 

IF(KF.EQ.l)  GO  TO  220 

IF( IDER(3).EQ.0>  GO  TO  240 

KA=1 

GO  TO  90 

CALCULATES  A   FOR  TRIMITE 

DO  230  IA=1,N1 

DIV(15,IA)=.5*(XI-X(IAH*{XI-X(IAII 

DIV(13,IA)=1.+DIV(15,IA)*DIV(10,IAI-3.0000000*(XI-X(IA|I 
It  IAI 


DRV2630I 
DRV2  640 
DRV2650 
DRV2660I 
DRV2670 
DRV2680I 
DRV2690I 
DRV2700I 
DRV2710J 
DRV2720I 
ORV27301 
DRV27401 
DRV2750t 
DRV2760( 
DRV2770( 
0RV2780( 
DRV2790( 
DRV2800( 
DRV2810I 
DRV2820( 
DRV2830( 
DRV2840( 
DRV2850( 
DRV2860( 
DRV2870C 
DRV2880( 
DRV2890( 
DRV2900C 
DRV2910C 
0RV2920( 
DRV2930C 
DRV2940C 
DRV2950C 
DRV2960( 
DRV2970C 
DRV2980C 
DRV2990(; 
DRV3000C 
DRV3010C 
DRV3020C 
ORV3030C 
DRV3040C 
DRV3050C 
0RV3060C 
DRV3070C 
DRV30800 
DRV3090C 
DRV31000 
DRV3U0G 
DRV31200 
ORV31300 
DRV31400 
DRV31500 
♦DIVC30RV31600 
DRV31700 
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DIV(14,IA|»(XI      -XUAII       ♦DIV(15,IA)*DIV(11.IAI 

IF{KF.EQ.1>    GO    TO    230                      '                J    OIVlll.IAl  DRV31800 

DIV(6,1)«DIV(1,IA|*DIV(1,IA|  DRV31900 

C  DRV32000 

C           CALCULATES    THIRD    DERIVATIVE    FOR   TRIMITE  0RV32100 

V*l 680.0000*0 1  V(  1,1  A)*DIV(  6, 11  DRV32300 

VA=15120.000*DIV(4,IAI*DIV(6,1|  DRV32400 
JB.226e0.0OO.DIVUfU,,D,vufIM,OIV(4fI„^^ 

DO    225    KFL    =1,3  DRV33100 

2"      ^'^•<F'-»^<U.1.KFL„60480.000      .HSQC  ^"foO 

M^!sUn!i!;nii*",i,,*,*cw,M,x,i,*w2',*,,Mw,»''»'»^-«».»»SwS»M 

230    CONTINUE  DRV33600 

IF(KF.EQ.l)  GO  TO  450  DRV33700 

IF  (  KB  .EQ.l)  GO  TO  240  DRV33800 

DVY(3)*DVY<3)/60480.000      *HSOC  DRV33900 

C                                      iWU  DRV34000 

C     CALCULATES  FOURTH  DERIVATIVE  FOR  TRIMITE  0^34200 

240    IF(IDER(4).EQ.O.  RETURN  DRV34300 

DO  251  IA=1,N1  DRV34400 

V=DIV(1,IA»«DIV(1,IA»  DRV34500 

VA  =  16  800. 000*V*DIV(1,IA)  DRV  34600 

VB=75600.000+V*DIV(4,IA»  DRV34700 

C0EF(IA,2,1»=VA*DIV(7,IAI*VB*DIV(10,IAI  nowl*!™ 

C0EF(IA.2.2I-VA*DIVC8,IAI*VB«0IV  11   A  SSX?t«?2 

C0EF<IA,2,3>  =  VA*DIV(9,IA)*VB*DIV  I    A  S!X«?S2 

IF(KB.EO.OJ  GO  TO  250               '   '  DRV35100 

DO  245  KFL=1,3  DRV35200 

245  j°E;si;j;^L«w«".2tMLi/isi»».oo  *hsqc  ski*™ 

251     CONTINUE  DRV35700 

IF  (  KB  .EQ.1»  RETURN  DRV35800 

DVYU)=DVY(4,/151200.00        *HSQC  DRV35900 

RETURN  DRV36000 

260   IF(MOLD  .GE.  MO  .AND.  KG  .EQ.l.  GO  TO  520  DRV362^ 

Nl«*  DRV36300 

GO  TO  30  DRV36400 

270   KA=1  DRV36500 

GO    TO    (290, 90, 1201, MO  DRV36600 

C  DRV36700 

C     CALCULATES  LUI  DRV36800 

C  DRV36900 

280   IF(KA.E0.2|  GO  TO  270  DRV37000 

290   DO  300  IA=1,N1  DRV37100 

DRV37200 
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300 
310 

320 

330 
3*0 

350 
360 

370 

380 
390 

400 


C 

c 
c 


DIV(16,IA)»ALPHA(1, I  A , 1) * Al PHA{ 1 , I A , 2 )* ALPHA(  1,  I A , 3  I /DI V( 6, I A ) 
GO  TO  (310, 380, 200), MOLD 
GO  TO  (360,340,320) ,M0 

INTERPOLATES  THE  SECONO  DERIVATIVE  USING  LAGRANGIAN 

GAM(3)«0. 

DO  330   IA*1,N1 

COEF( IA,1,3)=DIV(4,IA) 

IF(K8  .EQ.  1)  GO  TO  330 

GAM(3)  =  GAM(3H-DIV(4,IA)*FN(l,U) 

CONTINUE 

INTERPOLATES  THE  FIRST  DERIVATIVE  USING  LAGRANGIAN 

GAM  (2  )  =  0. 

DO  350   IA=1,N1 

COEF(IA,l,2)=DIV(2,IA I 

IF(KB  .EQ.  1)  GO  TO  350 

GAM(2)»       GAM(2)       +D I V(2 , I A)*FN( 1, IA ) 

CONTINUE 

INTERPOLATES  THE  FUNCTIONAL  VALUE  USING  LAGRANGIAN 

GAM(l)  *  0. 

DO  370   IA=1,N1 

COEF( IA,1,1)=DIV( 16, IA) 

IF(  KB  .EQ.  1  )  GO  TO  370 

GAM(l)  =  GAM(l)  ♦  DIV(16,IA)  *  FN(1,IA) 

CONTINUE 

RETURN 

DEVELOPS  COEFFICIENTS 


*OIV(  7,IA) 


DO  390   IA=1,N1 

DIV(7  ,IA)=-2.»DIV(3,IA) 

DIV(13,  IA)  =  1.+(XI-X(IA)  ) 

DIV(8,IA)=1. 

DIV(14,IA)=XI-X(IAI 

GO  TO  (43C,  420,400,  MO 


INTERPOLATES  SECOND  DERIVATIVE  USING  HERMITE 


410 


GAM 
DO 
V  =  2 
VA  = 
COE 
COE 
IF 
GAM 
1  IA) 
CON 


(3)  =  0. 

410  I 
•*(DIV 
4.*DIV 
F( IA,1 
F( IA,2 
(KG  .E 
(3)=GA 

TINUE 


A=1,N1 

(2,IA>*DIV(2,IA)+DIV(16,IA)*DIV(4,IA>) 

(16,IA)*DIV(2,  IA) 

,3)=V*DIV(13,IA)+VA*DIV(7,IA> 

,3)  =  V*DIV(  14,1  AKVA 

Q.  1)  GO  TO  410 

M(3)  +C0EF(IA,i,3)*FN(l,IA)+C0EF(IA,2,3)*FN( 


INTERPOLATES  FIRST   DERIVATIVE  USING  HERMITE 


DRV3730( 
DRV3740( 
DRV3750( 
0RV3760(j 
DRV3770C 
DRV3780d 
DRV379O0 
DRV38000 
DRV3810(J 
DRV38200 
ORV38300 
DRV38400 
DRV38500 
ORV3860C 
DRV  38700 
DRV38800 
DRV38900 
DRV39000 
DRV39100 
DRV39200 
DRV39300 
DRV39400 
DRV39500 
DRV39600 
DRV39700 
DRV39830 
DRV39900 
DRV40000 
DRV40100 
DRV4020O 
DRV40300 
DRV40400 
DRV40500 
DRV40600| 
DRV40700 
DRV40800;! 
DRV4090Q 
DRV41000, 
DRV41100 
DRV41200 
DRV41300 
DRV41400 
DRV41500 
DRV416Q0 
DRV41700 
DRV41800 
DRV41900 
DRV42000 
DRV42100 
DRV42200 
2,DRV42300 
DRV42400 
DRV42500 
DRV42600 
DRV42700 


99 


420      GAM(2I«0. 

DO    425       IA»l,Nl 

V*2.*DIV(16,IAI*DIV<2,IAI 

VA«DIV(16,IAI»DIV(16,IA» 

C0EF(IA,1,2»*V*DIV(13,IAH-VA*DIV<7.IAI 
C0EF(IA,2,2»xV*DIV<l4,IA»+VA  l',IAI 

IF(KB    .  EQ.    1)    GO    TO   425 

CONTINUE    GAM<2'    *   C0EF(IA»1»2>    *FN(l,IA) 


425 


♦    C0EF(IA,2,2»*FN(2,IA| 


43  0 


440 


450 


C 
C 
C 


INTERPOLATES    THE    FUNCTIONAL   VALUE    USING   HERMITE 

GAM(1J=0. 

DO    440       IA"1,N1 

VA=DIV(16,IA)*DIV(16,IA» 

COEFt IA,i,l|=VA*DIV(13, IA  I 

COEF( IA,2,1 J=VA       *DIV(14,IAI 

IF(K8    .  EQ.    1  I    GO    TO   44C 

CONTINUE    GAM<U    *    C0EF(IA»1*1'*FN(l,IAI 

RETURN 

GO    TO       (500, 480,460), MO 


COEFUA,2,l>*FN(2,IA) 


460 


470 


INTERPOLATES    SECOND    DERIVATIVE    USING    TRIMITE 

GAM(3)=0. 

DO    470       IA*1,N1 

VB=DIV(16,IAJ*DIV(16,IAI 
VA=6.  0000000* VB*D IV (2, 1  A) 

ve=^im:TA?,V,4'U,t6-Coooaoc*0,v(l6'IA,'OIV'^'*»*°'v,2..A, 
cnlt   "•!!•?!"'•<>""  ".'A. ♦vA.ome     A MHo  V  11     A 

S«8"M:n,"S0Tol"iu,w"wy",,u,*vwoM,»:»» 

CONTINUE 


♦    C0EF(IA,2,3)*FN(2, IAI 


480 


490 


INTERPOLATES  FIRST   DERIVATIVE  USING  TRIMITE 

GAM(2>=0. 

DO  490   IA=1,N1 

VB=DIV(16,IAI*DIV(16,IA» 

V=3. 0000000* VB*DI V(2,IAI 

VA=0IV(16,IA)*VB 

COEF(  IA,1,2I  =  V»DIV(13,IAH-VA*DIV(7,IAI 

C0EF( IA,2,2,  =  V*DlV(l4,IAH-VA*DIV<8,IA) 

rc^iIA:3,2,=V¥DIV(15*IA'+VA*0^(9,IA) 
IF(KB  .EQ.  1 )    GO  TO  490 

GAM(2>  =  GAM(2I  ♦  C OEF ( I  A, 1 , 2 )*FN( 1 , I  A) 
I  ♦  C0EF(IA,3,2)*FN(3,IAI 
CONTINUE 


COEF(IA,2,2l*FN(2,IAI 


DRV42800 

DRV42900 

DRV43000 

DRV43100 

DRV43200 

DRV43300 

DRV43400 

DRV43500 

ORV43600 

DRV43700 

DRV43800 

DRV43900 

DRV44000 

0RV44100 

DRV44200 

DRV44300 

DRV44400 

DRV44500 

DRV44600 

DRV44700 

DRV44800 

DRV44900 

DRV45000 

DRV45100 

DRV45200 

DRV45300 

0RV45400 

DRV45500 

DRV45600 

DRV45700 

DRV45800 

0RV45900 

0RV46000 

DRV46100 

DRV46200 

DRV46300 

DRV46400 

DRV46500 

DRV46600 

DRV46700 

DRV46800 

DRV46900 

DRV47000 

DRV47100 

DRV47200 

DRV47300 

DRV47400 

DRV47500 

DRV47600 

DRV47700 

DRV47800 

DRV47900 

DRV48000 

0RV48100 

DRV48200 


l.OQ 


INTERPOLATES    THE    FUNCTIONAL   VALUE    USING    TRIMITE 


5CO      GAM(1I»0. 

DO    510       IA»1,N1 

V=0IV(16,IA>*DIV(16,IA>*DIV(16,IA) 

COEF{ IA,l,ll=V*DIV(13,IAI 

COEF( IA,2,1I  =  V*DIV(  14.IAI 

COEF{ IA,3tl )=V*DIV(15,IA» 

IFCKB    .EQ.    1)    GO   TO    510 

GAM(l)    =    GAM(l)    +    COEF(IA,i,l)*FN(l,IAI 
1    ♦    COEF< IAt3tll*FN(3t IAI 
510      CONTINUE 

RETURN 
520      WRITE(6,1C?01) 
10001    FORMAT!*     INVAL ID   ORDERS • I 

RETURN 

END 


♦  C0EF(IA,2,1H"FN<2,  I  A) 


DRV4830C 
DRV4840C 
DRV4850C 
DRV4860C 
DRV4870C 
DRV4880C 
DRV4890C 
DRV4900C 
DRV4910C 
DRV^920C 
DRV4930C 
DRV49400 
DRV4950C 
DRV49600 
DRV49700 
DRV49800 
DRV49900 


C 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


SUBROUTINE  DECOMP <N,UL, J5 ,K7 ) 
IMPLICIT  REAL  +  8(A-H,0-ZI 

THIS  SUBPROGRAM  DECOMPOSES  THE  MATRIX  INTO  UPPER  AND  LOWER 
TRIANGULAR  MATRICES.   THESE  Ha D  THE  OPERATIONS  PERFORMEO  IN 
THE  FIRST  HALF  OF  A  GAUSSIAN  ELIMINATION  PROCEDURE, 

THE  VARIABLES  IN  THIS  SUBPROGRAM  ARE 
UL  THE  ARRAY 

N  THE  NUMBER  OF  VARIABLES 

J5  A  RETURN  CODE  1  IF  SUCCESSFUL,   -1  IF  NOT 

K7  THE  NUMBER  OF  DIAGONALS  IN  THE  MATRIX 

K3  POINTS  TO  POSITION  OF  THE  MAIN  DIAGONAL 

K2  NUMBER  OF  SUBDIAGONALS 

THE  REST  OF  THE  VARIABLES  ACT  AS  COUNTERS  OR  AS  POINTERS 

TO  POSITIONS  IN  UL 

DIMENSION  UL(K7,1» 

COMMON  /CNT/  ISTEP, NUMDE, NUMSO 

NUMDE=NUMDE+1 

J5=l 

K2=K7/2 

K3=K2+l 

DO  10  I=K3,N 

K  =  I  ♦  1 

I1=I-K2 

DO  10  J=1,K2 

K  =  K  -  1 

IF(DABS(ULtK3,Il»I.LT,i.D-2C*)   GO  TO  15 

CALCULATES  THE  FACTOR  USED  TO  SUBTRACT  ONE  EQUATION  FROM 

ANOTHER 


DCP10000 

DCP10100 

DCP1D200 

DCP10300 

DCP10400 

DCP10500 

OCP10600 

DCP10700 

DCP10800 

DCP10900 

DCP11000 

DCP11100 

DCP11200 

DCP11300 

DCPlUOOi 

DCPU500 

DCP1160C 

DCP11700 

DCP11800 

DCP11900 

DCP12000 

DCP12100 

DCP12200 

DCP12300 

DCP12400 

DCP12500 

DCP12600 

DCP12700 

DCP12800 

DCP12900 

DCP13000 

DCP13100 

DCP13200 


101 


c 
c 
c 


10 


c 
c 
c 
c 


20 
15 


UL(J,K»»UL(J,KI/Ul(K3,m 
J1»J*1 

J2»J*   K2 

LI  -  K3 

DO  13  L  «  JitJ2 

LI  -  LI  ♦  1 

SUBTRACTS  CORRESPONDING  ELEMENTS  FROM  THE  TWO  EQUATIONS 

UL(L,K>*UL(L,K)-UL(L1,UI*JLCJ,KI 
IF(K2.EQ.l)  RETURN 

DO  20  I»2,K2 

11  »  N  ♦  I 
K5  »  U  -  K3 

12  *  K2  -  I 
DO  20  J«I,K2 
K  *  U  -  J 

IF(DABS(UL(K3,K5M.LT.l.D-20»   GO  TO  15 

UL(J,K>=UL(JfK)/UL(K3,K5) 

J1»J+1 

J2  »  12  ♦  Jl 

LI  *  K3 

DO  20  L»J1,J2 

LI  »  LI  ♦  1 

UL(L,K»=UL(L,K»-UL(L1,K5I*UL(J,KI 

J5=-l 

RETURN 

END 


DCP13300 

DCP13400 

DCP13500 

DCP13600 

DCP13700 

DCP13800 

DCP13900 

DCP14000 

DCP14100 

DCP14200 

DCP14300 

DCP14400 

DCP14500 

0CP146C0 

DCP14700 

DCP14800 

DCP14900 

DCP15000 

DCP15100 

DCP15200 

DCP15300 

DCP15400 

DCP15500 

DCP15600 

OCP15700 

DCP15800 

DCP15900 

DCP16000 

DCP16100 

DCP16200 

DCP16300 

DCP16400 

DCP16500 


SUBROUTINE  SOLVEJ N.UL, B, X, K7,  ID,K> 
IMPLICIT  REAL  *8(A-H,0-ZI 

THIS  PROGRAM  USES  THE  DECOMPOSED  MATRIX  FROM  DECQMP  Tn 
PERFORM  THE  GAUSSIAN  ELIMINATION  T° 


THE 


VARIABLES 
UL 

N 
K7 
K3 
K2 

B 

X 
ID 

K 


IN  THIS  SUBPROGRAM  ARE 
THE  ARRAY 

THE  NUMBER  OF  VARIABLES 
THE  NUMBER  OF  DIAGONALS  IN  THE  MATRIX 

iS  S  £5?JIKsTHE  MA,N  0UG0N4L 

thI  "suLTNJ6cVrm°R  °F  ™E  E0UAr,0N  I0  BE  S0L'E0 

THE  NUMBER  OF  ROWS  IN  THE  RESULT  VECTOR 
:TuL??W  IN  RESULT  VECT0R  WHERE  ANSWER  IS  PLACED 


THE  REST  OF  THE  VARIABLES  ACT  a"s  COUNTERS  OR  A  S  PINTERS 


SLV10000 

SLV10100 

SLV10200 

SLV10300 

SLV10400 

SLV10500 

SLV10600 

SLV10700 

SLV10800 

SLV10900 

SLV11000 

SLV11100 

SLV11200 

SLV11300 

SLV11400 

SLV11500 

SLV11600 


102.  . 


TO  POSITIONS  IN  UL 

DIMENSION    UL(K7,i),X< ID, II, Bill 
COMMON    /CNT/    ISTEP,NUMDE,NUMSO 
NN    =    N    ♦    1 
NUMSO*NUMSO*l 

FIRST    THE    RESULT    VECTOR    IS    INITILIZED    TO    THE    CONSTANT    VECTOR 

D01I=*1,N 
X(K,I »=B(  I) 


K2=K7/2 
K3=K2+1 
K4=K2-1 
IF(K2.EQ.  1>  G0T09 


THIS  SECTION  PERFORMS  THE  OPERATIONS  OF  THE  FIRST  PART  OF  THE 
GAUSSIAN  ELIMINATION  ON  THE  CONSTANT  VECTOR 

D02  I  =  2,K2 
Il»  1-1 
J3  =  K3  -  I 
DO  2  J=1,I1 
J2  =  J3  ♦  J 
2    X(K,I >=X(K,I)-X(K,J)*UL(J2  ,11 

9  DO  3  I=K3,N 
D03  J  =1,K2 
J2=I-J 
J3=K3-J 

)     X(K,I >=X(K, I >-X<K,J2)*UL(J3,I I 
IF(K2.GT.l)  GO  TO  W 
X(K,N)=X(K,NI/UL(2,N» 
GO  TO  5 

4  X(K,N)=X(K,NI/UL(K3,N» 

THIS  SECTION  HANDLES  THE  BACK  SUBSTITUTION 

DO  6  J  =  1,K* 

NJ=N-J 

DO  7  1=1, J 

KJ=K3+I 

NJ1=NJ+I 

7  X(K,NJ>=X(K,NJ>-UL(KJ,NJ)*X(K,NJ1I 
6    X(K,NJ)  =X(K,NJI/UL(K3,NJ» 

5  DO  10  IBK=K3,N 
I   =  NN  -  IBK 
DO  8  J  =  1,K2 
JK=K3+J 

IJ=  I+J 

8  X{K,n  =  X(K,I)-UL(  JK,  I»*X(K,  IJ) 

10  X(K,I  >=X(K,I)/UL(K3,II 
RETURN 

END 


! 


SLV1170Q 

SLV1180( 

SLV1190I 

SLV1200< 

SLV1210( 

SLV12200; 

SLV12300 

SLV12400 

SLV12500 

SLV12600 

SLV12700 

SLV12800! 

SLV12900 

SLV13000 

SLV13100 

SLV13200 

SLV13300 

SLV13400 

SLV13500 

SLV13600 

SLV13700 

SLV13800 

SLV13900 

SLV14000 

SLV14100 

SLV14200 

SLV14300 

SLV14400 

SLV14500 

SLV14600! 

SLV14700 

SLV14800 

SLV14900 

SLV15000 

SLV15100 

SLV15200! 

SLV15300 

SLV15400' 

SLV15500 

SLV15600 

SLV15700 

SLV15800 

SLV15900 

SLV16000 

SLV16100 

SLV16200 

SLV16300 

SLV16400 

SLV16500 

SLV16600 

SLV16700 

SLV16800 

SLV16900 

SLV17000 

SLV17100 


103 


C 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


2. 
3. 


SUBROUTINE  OIFFUN( X,Y,T,D IVER.MO, INT ,N, IDER,ERR,OY,L.  12, L3 I 

THE  PURPOSE  OF  THIS  PROGRAM  IS  TO  PROVIDE 

1.  THE  TIME  DERIVATIVE  AS  THE  RESULT  OF  THE 

DIFFERENTIAL  EQUATION  DURING  INTEGRATION 

THE  TIME  DERIVATIVE  WHEN  STARTING 

THE  TIME  DERIVATIVE  FOR  THE  CALCULATION  OF  THE 

JACOB  I  AN 
IN  THE  FIRST  CASE  L  *  L2  «  0 

IN  THE  SECOND  CASE  L*1,L2»0  AND  THE  PROGRAM  ESTIMATES  THE  TIMF 
POINTsIIVE  °F  ™E  B0UN°ARY  P0INTS  8V  EXTRAPOLATION  FROM  INTERIOR 
ARETBYpIsSED.CASE  LC°fl2'1    AN°  ™E  ERR0R  ESTIMATION  CALCULATIONS 
THE  VARIABLES  IN  THIS  SUBPROGRAM  ARE 


X 

Y 

T 
DIVER 

MO 
INT 

N 

IDER 

ERR 

DY 

L 

L2 

DV 

DVP 

THE  REST 


HOLDS  THE 
HOLDS  THE 
TIME 

CONTAINS 
INTERPOLA 
ORDER  OF 
PROVIDE  P 
STARTING 
NUMBER  OF 
CONTAINS 
CONTAINS 
CONTAINS 
INDICATES 
INDICATES 
TEMPORARI 
TEMPORARI 
OF  THE  VARIABLE 


PARTITION  POINTS 
FUNCTIONAL  VALUES 


THE  COEFFICIENTS  OF  ERROR  ESTIMATE 

TING  POLYNOMIALS 

THE  METHOD 

OINTS    OF    ANALOG    f INTC 1*1 II    AND    THE 

POINT  OF  ERROR  ESTIMATE  ANALOG  (INT(2tIII 

EQUATIONS 
NAMES  OF  DERIVATIVES  FOR  ALL  ORDERS 
ERROR  ESTIMATE  VECTOR 
THE  TIME  DERIVATIVE 

IF  FINDING  TIME  DERIVATIVE  AT  BOUNDARY 

JACOBIAN    CALCULATION 
LY    STORES    SPACE    DERIVATIVES 
LY    STORES    ERROR    ESTIMATE    FOR    DV 
S    ACT    AS    COUNTERS    OR    POINTERS. 


IMPLICIT  REAL  *  8  (A-H,Q-Z» 
DIMENSION  SV(6» 

1o^^«:!Um:!:o?;::^j;^:i:,dwer,9,mo'i,',nt,2'i,'ider«4'- 

DIMENSION  XT(5),XS(5I,YT<3,5I,YS<3,5I 

COMMON  /DVR/  DVY(  4  »  ,COEF  (  4,  3  f  31  ,  ALPHA!  2,4,  3  >,  DI  V(  16,  5  ) 

THIS  SECTION  DETERMINES  THE  TIME  DERIVATIVES  AT  THE  BOUNDARIES 

IF(L3.EQ.0I  GO  TO  4 
DO  3  IA*l,MO 
SV(IAJ*Y(1,IA,1I 

SV(IA*3I=Y(1,IA,NI 

CONTINUE 

CALL  BDRY(XfY,T,MO,N,DY) 

DO  5  IA  =  l,MO 


DFNIOOOO 

DFN10100 

DFN10200 

DFN10300 

DFN10400 

DFN10500 

DFN10600 

DFN10700 

DFN1Q800 

DFN10900 

DFNllOOO 

DFNlllOO 

DFN11200 

DFN1130O 

DFN11400 

DFN11500 

DFN11600 

DFN11700 

DFN11800 

DFN11900 

DFN12000 

DFN12100 

DFN1220O 

DFN12300 

DFN12400 

DFN12500 

DFN12600 

DFN12700 

DFN12800 

DFN1290C 

DFN13000 

DFN1310O 

DFN13200 

DFN13300 

DFN13400 

DFN13500 

DFN13600 

DFN13700 

DFN13800 

DFN13900 

DFN14000 

DFN14100 

DFN14200 

DFN143C0 

DFN14400 

DFN14500 

DFN14600 

DFN14700 

DFN14800 

DFN14900 


ioU 


c 
c 
c 


10 


20 


C 

C 

c 
c 


25 


30 


40 


50 


60 


ERR(IA,1  I  »  DY(  IA,1  I 
ERR(IA,N>  »  DY(  IA,NI 
CONTINUE 

THIS    SECTION    DETERMINES    WHICH   POINTS    ARE    USED    IN    THE    DERIVATIVE 

NM1*N-1 

DO       100       IA    *    2.NM1 

IG»IA 

IF(  INTd,  IAI.LT.O)    GO   TO    10 

IB=IA+INT(1,IA» 

IC=IA-1 

DETERMINES  STARTING  POINT  FOR  ERR  AND  THE  POWERS  OF  SPACING 

GO  TO  20 

IB=IA+1 

IC=IA*INT(1,IAI 

ID=INT(2,IAI 

XI=X(IAI-X{ IC) 

IF  (XI    .  EQ.    XITMP      .AND,       MO    .  EQ.    MOLD!    GO    TO    25 

MOLD    =    MO 

XITMP    *    XI 

XIS=XI*XI 

IF    (MO    .LT.    2)    GO    TO    25 

XIC=XI*XIS 

IF    (MO    .LT.    31    GO    TO   25 

XIF=XI*XIC 

GO    TO    (30,70,80  I, MO 

SETS  PARAMETERS  AND  FORMS  DERIVATIVES  FOR  LAGRANGIAN 

IF(IDER(1>.EQ.1I    DV(1I=(Y(1,1,IBI-Y(1,1,IC)I/(2.*XII 
IFC I0ERC2l.EQ.il    DV( 2 1  = ( Y( 1 , 1 ,1 B 1-2. *Y( 1, 1 , IA >+Y( 1 , 1, IC ) I /XIS 
DVP(1I=0. 
DVP(2I=0. 

USES    COEFFICIENTS    OF    INTERPOLATING    POLYNOMIAL    AND    CALCULATES 
ERROR    DERIVATIVES. 

IF    (L2    .EQ.    1)    GO    TO   65 
DO    50    IE=1,M0 
DO    5D    IF=1,4 
IH=ID+IF-1 

IFl  IDER(M0|.EQ.1J    DVP( MOI =D I VER ( I F, I E , I A  I *Y( 1 , IE, IH )+DVP( M3 t 
IF(IDER(MO+ll.EQ.il    DVP ( M0*1I =D IVER ( I F+4, IE, I Al *Y( 1 , I E, IH l+DVP( MO 
1+1  I 
CONTINUE 

IF(MO.NE.l)    GO    TO    60 
IFC  IDERC2l.EQ.il    DVP(2I=D  IVER  (9, 1    ,  I  A)*Y(  1, 1,  ID+4  KDVPl  21 

USES    DIFFERENTIAL    EQUATION   TO    OBTAIN    TIME    DERIVATIVES    AND   ERRORS 

IF(  IDER(M0I.EQ.1I  DVP( MOI =DVP(MO I *DV( MO  I 

IF(  IDER(M0*1I.EQ. II  DVP  (  MO*  1  I  =DVP  (MOU  I  +DVC  MO+1  I 


DFN15000 

OFN15100 

DFN15200 

DFN15330 

DFN15400 

DFN15500 

DFN15600 

DFN15700 

DFN15800 

DFN15900 

DFN16000 

DFN16100 

DFN16200 

DFN163Q0 

DFN16400 

DFN16500 

DFN16600 

DFN16700 

DFN16800 

DFN16900 

DFN17000 

0FN17100 

DFN17200 

DFN17300 

DFN17400 

DFN17500 

DFN17600 

DFN17700 

DFN17800 

DFN17900 

DFN18000 

DFN18100 

DFN18200 

DFN18300 

DFN18400 

DFN18500 

DFN18600 

DFN18700 

DFN18800 

DFN18900 

DFN19000 

DFN19100 

DFN19200 

DFN19300 

DFN19400 

DFN19500 

DFN19600 

DFN19700 

DFN198D0 

DFN19900 

DFN20000 

DFN20100 

DFN20200 

DFN20300 

DFN20400 


105 


65 


100 


103 
105 


110 


120 


CALL  UFUN(M0,XUA),Y<lfl,IA»fT,DVP,ERRU,IA»,5l 

CALL  UFUN(MO,X(IAJ,Y(lfltIA),T,DV,DYIl,IA|,5l 

DO  100  IB»1,M0 

IF(L2.EQ. II  GO  TO  100 

ERRflBtlAI  -  DY(IB,IA|  -  ERRUB,IAI 

CONTINUE 

IFIL3  .  EQ.  0)  GO  TO  135 

DO  103  IA«1,M0 

Y(l,IAtl)*SV(  IA  I 

Y(1,IA,N   l»SVUA*3> 

CONTINUE 

0  )  RETURN 


IF 
XI 
XN 
DO 
II 
12 


C 


L  .EQ. 

•  xui 

'  X(N) 

110  I  *   1,5 
'    I  ♦  1 
»  N  -  I 
XT(II  »  X(Il» 

xsm  *  xii2» 

00   110  J  *  1,M0 

YT(J,II  »  DYU.Il) 

YS( J,I)  «  DY( J, 12) 

CONTINUE 

CALL  DERIVUT.YT,  I0ER 

CALL  DERIV(XS,YS, IDER 

DO   120  I  ■  1,M0 

DY(  I» II  *  DC  II 

DY(  I,NI  *  ER(I) 

CONTINUE 

RETURN 


,M0, 0, 1,D,M0,X1,D,  XII 
»M0,0, l,0,MOtXN,ER,XII 


SETS  PARAMETERS  AND  FORMS  DERIVATIVES  FOR  HERMITE 

70    DV(1I=Y(1,2,IGI 
DVP(1I*DV(1I 

IF(  IDER(2).EQ.CI  GO  TO  75 
DVP(2I=0. 

ic^i/(2!*xn(1,1'IB,"2'*Y{1,1,IG,4Y(1,1,ICn/,(IS',Y,1,2,IB,"Y(U2' 

75   IF( IDER(3).EQ.0I  GO  TO  UO 
DVP(3»>=0. 

DV(3l=7.5000C0C*(Y(ltl,IBI-Y(lfl,IC)l/XIC-1.500000D*(Y(l,2.IBI 
1  ♦8.*Y(1,2,IG)  ♦  Y(1,2,ICII/XIS 

GO  TO  40 

SETS  PARAMETERS  AND  FORMS  DERIVATIVES  FOR  TRIMITE 

80    DVIll=Y(l,2tIG» 
DVP(1)*DV(1J 
DV(2I=Y(1,3,IG» 
DVP(2I=DV(2I 

IF( IDER(3). EQ.OI  GO  TO  85 
DVP(3I»0. 

1  +Y(lt2, IC»l*18.000000»Y(lf2,IGII/XIS+.37500000*(Y(l,3,IB»  - 


OFN205OO 
DFN20600 
DFN20700 
DFN20800 
DFN20900 
DFN21000 
DFN21100 
DFN21200 
DFN21300 
DFN21400 
DFN21500 
0FN21600 
DFN21700 
DFN21800 
DFN21900 
DFN22000 
DFN22100 
DFN22200 
OFN22300 
DFN22400 
DFN22500 
DFN22  600 
DFN22700 
DFN22800 
DFN22900 
0FN23000 
DFN23100 
DFN23200 
DFN23300 
DFN23400 
DFN23500 
DFN23600 
DFN23700 
DFN23800 
DFN23900 
DFN24000 
DFN24100 
IDFN24200 
0FN24300 
DFN24400 
DFN24500 
DFN24600 
DFN24700 
DFN24800 
DFN24900 
DFN25000 
0FN25100 
DFN25200 
DFN25300 
DFN25400 
DFN25500 
DFN25600 
DFN25700 
OFN25800 
DFN25900 


106 


c 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


2   Y(lf3tICII/XI  DFN26O00 

85   IF{  IDERm.EQ.O)  GO  TO  40  DFN26100 

DVP(4)«0.  DFN26200 

DV(4)«72.000000*(Y(1,1,IBI-2.*Y(1,1,IG)«-Y(1,1,IC)  )/XI F-19. 500000*  DFN26300 

l(Y{l,2t IBJ-Y(l,2,IC)>/XIC*1.5COOOOa*(Y(l,3,IB)-24.00000C*Y(l,3,IG)DFN26400 

2  ♦Y(1,3,ICM/XIS  DFN26500 

GO  TO  40  DFN26600 

END  DFN26700 


SUBROUTINE    01 FSUB (N ,T ,Y .SAVE, H, HMIN, HMAX, EPSt YMAX, ERROR ,KFLAG , 
USTART,MAXDER,PSAVE,K7,HMY,X,DIVER,INT,IDER,ERR    , MO , I OOUB, I QQQC, 
2HSAV,PEQ,ER0RA,ITPSC) 

IMPLICIT  REAL  *  8(A-H,P-Z» 

REAL    *    4    PEPSH,PRl,PR2tPR3,PERTST 

THE  PARAMETERS  TO  THE  SUBROUTINE  DIFSUB  HAVE 
THE  FOLLOWING  MEANINGS*. 


N 


SAVE 


HMIN 


HMAX 


THE  NUMBER  OF  FIRST  ORDER  DIFFERENTIAL  EQUATIONS.   N 
MAY  BE  DECREASED  ON  LATER  CALLS  IF  THE  NUMBER  OF 
ACTIVE  EQUATIONS  REDUCES.  BUT  IT  MUST  NOT  BE 
INCREASED  WITHOUT  CALLING  WITH  JSTART  »  0. 

THE  INDEPENDENT  VARIABLE. 

AN  8  BY  N  ARRAY  CONTAINING  THE  DEPENDENT  VARIABLES  AND 
THEIR  SCALED  DERIVATIVES.   Y(J+1,II  CONTAINS 
THE  J-TH  DERIVATIVE  OF  Y(I)  SCALED  BY 
H**J/FACTORIAL(J)  WHERE  H  IS  THE  CURRENT 
STEP  SIZE.  ONLY  Yd,  II  NEED  BE  PROVIDED  BY 
THE  CALLING  PROGRAM  ON  THE  FIRST  ENTRY. 

IF  IT  IS  DESIRED  TO  INTERPOLATE  TO  NON  MESH  POINTS 
THESE  VALUES  CAN  BE  USED.   IF  THE  CURRENT  STEP  SIZE 
IS  H  AND  THE  VALUE  AT  T  +  E  IS  NEEDED,  FORM 
S  =  E/H,  AND  THEN  COMPUTE 
NQ 
YUMT+EI  »   SUM   YU  +  1,1  )*S**J 
J=0 

A  BLOCK  OF  AT  LEAST  12*N  FLOATING  POINT  LOCATIONS 
USED  BY  THE  SUBROUTINES. 

THE  STEP  SIZE  TO  BE  ATTEMPTED  ON  THE  NEXT  STEP. 
H  MAY  BE  ADJUSTED  UP  OR  DOWN  BY  THE  PROGRAM 
IN  ORDER  TO  ACHEIVE  AN  ECONOMICAL  INTEGRATION. 
HOWEVER,  IF  THE  H  PROVIDED  BY  THE  USER  DOES 
NOT  CAUSE  A  LARGER  ERROR  THAN  REQUESTED,  IT 
WILL  BE  USED.   TO  SAVE  COMPUTER  TIME,  THE  USER  IS 
ADVISED  TO  USE  A  FAIRLY  SMALL  STEP  FOR  THE  FIRST 
CALL.   IT  WILL  BE  AUTOMATICALLY  INCREASED  LATER. 

THE  MINIMUM  STEP  SIZE  THAT  WILL  BE  USED  FOR  THE 
INTEGRATION.   NOTE  THAT  ON  STARTING  THIS  MUST  BE 
MUCH  SMALLER  THAN  THE  AVERAGE  H  EXPECTED  SINCE 
A  FIRST  ORDER  METHOD  IS  USED  INITIALLY. 

THE  MAXIMUM  SIZE  TO  WHICH  THE  STEP  WILL  BE  INCREASED 


DSB10000 

DSB10100 

DSB10200 

DSB10300 

DSB10400 

DSB10500 

0SB10600 

DSB 10700 

DSB10800 

DSB10900 

DSB11000 

DSBIHOO 

DSB11200 

DSB11300 

DSB11400 

DSB11500 

DSB11600 

DSB11700 

DSB11800 

OSB11900 

DSB12000 

DSB12100 

0SBI2200 

DSB12300 

DSB12400 

DSB12500 

DSB12600 

DSB12700 

DSB12830 

DSB12900 

DSB13000 

DSB13100 

DSB13200 

DSB13300 

DSB13400 

0SB13500 

DSB13600 

DSB13700 

DSB13800 

DSB13900 

DSBKOOO 

DSB14100 
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EPS 


MF 


YMAX 

ERROR 
KFLAG 


JSTART 


MAXDER 


AN  ARRAY  OF 
OF  EACH  Y 


THE  ERROR  TEST  CONSTANT,   SINGLE  STEP  ERROR  FSTimatpc 
DIVIDED  BY   YMAXUI   MUST  BE  LESS  THAN *?*, ESTlMATES 

IwUSTEO^^^ErSllJ*  STEP  AND/°R  °RDER  IS 
THE  METHOD  INDICATOR.  THE  FOLLOWING  ARE  ALLOWED.. 

0  AN  ADAMS  PREDICTOR  CORRECTOR  IS  USED. 

1  A  MULTI-STEP  METHOD  SUITABLE  FOR  STIFF 

lYntTlrLlSr°SBD'     IT  WILL  ALS0  WORK  FOR 
NON  STIFF  SYSTEMS.   HOWEVER  THE  USER 

MUST  PROVIDE  A  SUBROUTINE  PEDERV  WHICH 

EVALUATES  THE  PARTIAL  DERIVATIVES  OF 

THE  DIFFERENTIAL  EQUATIONS  WITH  RESPECT 

TO  THE  Y«S.    THIS  IS  DONE  BY  CALL 

PEDERV{T,Y,PW,M>.  PW  IS  AN  N  BY  N  ARRAY 

WHICH  MUST  BE  SET  TO  THE  PARTIAL  OF 

THE  I-TH  EQUATION  WITH  RESPECT 

TO  THE  J  DEPENDENT  VARIABLE  IN  PW(I.J). 

PW  IS  ACTUALLY  STORED  IN  AN  M  BY  M 

ARRAY  WHERE  M  IS  THE  VALUE  OF  N  USED  ON 

THE  FIRST  CALL  TO  THIS  PROGRAM. 

2    THE  SAME  AS  CASE  1,  EXCEPT  THAT  THIS 

SUBROUTINE  COMPUTES  THE  PARTIAL 

DERIVATIVES  BY  NUMERICAL  DIFFERENCING 

OF  THE  DERIVATIVES.  HENCE  PEDERV  IS 

N  LOCATIONS  WHICH  CONTAINS  THE  MAXIMUM 
SEEN  SO  FAR.   IT  SHOULD  NORMALLY  BE  SET  TO 

ANON4"R^E?FE^0E1NM^IcSHHCH^^C^T.MNS  THE  """"» 
i  COMPLETION  CODE  WITH  THE  FOLLOWING  MEANINGS.. 
THE  STEP  WAS  SUCCESSFUL 
THE  STEP  WAS  TAKEN  WITH  H  =  HMIN,  BUT  THE 

REQUESTED  ERROR  WAS  NOT  ACHIEVED. 
THE  MAXIMUM  ORDER  SPECIFIED  WAS  FOUND  TO 

BE  TOO  LARGE. 
CORRECTOR  CONVERGENCE  COULD  NOT  BE 

ACHIEVED  FOR  H  .GT.  HMIN. 
THE  REQUESTED  ERROR  IS  SMALLER  THAN  CAN 
BE  HANDLED  FOR  THIS  PROBLEM. 
AN  INPUT  INDICATOR  WITH  THE  FOLLOWING  MEANINGS.. 
-1   REPEAT  THE  LAST  STEP  WITH  A  NEW  H 
C   PERFORM  THE  FIRST  STEP.   THE  FIRST  STEP 
MUST  BE  DONE  WITH  THIS  VALUE  OF  JSTART 
SO  THAT  THE  SUBROUTINE  CAN  INITIALIZE 
ITSELF. 

.c-r.oT  +1   TAKE  A  NEW  STEP  CONTINUING  FROM  THE  LAST. 
ir  TIS  ^\™  NQ'  THE  CURRENT  ORDER  OF  THE  METHOD 

AT  EXIT.   NQ  IS  ALSO  THE  ORDER  OF  THE  MAXIMUM 
DERIVATIVE  AVAILABLE.  "AXIMUM 

THE  MAXIMUM  DERIVATIVE  THAT  SHOULD  BE  USED  IN  THE 
METHOD.   SINCE  THE  ORDER  IS  EQUAL  TO  THE  HIGHEST 

SI  ltllTE"SE0'    THIS  RESTRICTS  THE  ORDER.  IT  mIsT 
BE  LESS  THAN  8  FOR  ADAMS  AND  7  FOR  STIFF  METHODS. 


♦  1 
-1 

-2 

-3 

-4 


DSB14200 

DSB143Q0 

DSB14400 

DSB14500 

DSB14600 

DSB14700 

DSB14800 

DSB14900 

DSB15000 

DSB15100 

DSB15200 

0SB15300 

DSB15400 

DSB15500 

DSB15600 

DSB15700 

DSB15800 

DSB15900 

DSB16000 

0SB16100 

DSB16200 

DSB16300 

0SB16400 

DSB16500 

DSB16600 

DSB16700 

0SB16800 

DSB16900 

DSB17000 

DSB17100 

DSB17200 

DSB17300 

DSB17400 

DSB17500 

DSB17600 

0SB17700 

DSB17800 

DSB17900 

DSB18000 

DSB18100 

DSB18200 

DSB18300 

0SB18400 

DSB18500 

DSB18600 

DSB18700 

DSBI8800 

DSB18900 

DSB19000 

DSB19100 

DSB19200 

DSB19300 

DSB19400 

DSB19500 

DSB19600 


PSAVE 

HMY 

HSAV 

PEG 
ERR 
ERORA 

COMMON 
DIMENS 
DIMENS 
DIMENS 
DIMENS 
I  A(8I, 
DIMENS 
DIMENS 
DATA  X 


A  BLOCK  OF  AT  LEAST  N**2  FLOATING  POINT  LOCATIONS, 

STORES  THE  ASYMPTOTIC  ERROR  VECTOR 

SAVES  THE  PREVIOUS  VALUES  OF  THE  ASYMPTOTIC  ERROR 

AND    IS    USED    TO    STORE    DERIVATIVES 
STORES    THE    JACOBIAN 

CONTAINS    THE    SPATIAL    ERROR    ESTIMATE 
CONTAINS    THE    ESTIMATED,       INTRODUCED    ERROR 

/DVR/    DVY(4) ,COEF(4,3,3),ALPHA(2,4,3l,DIV(16,5l 
ION    XT(5t2lvXI(2l 
ION    YT(3,5> 

ION    V    <1 I ,D 

ION    Y(5,15  t 

PERTST(7,3J 

ION    ERORAd 

IONHMY(2,l»  ,HSAV(  3,1) ,PEQIK7,1» 

ITMP    /O.Q/ 


IVERCDt  INTClIt  IDERdlfERRCll      ,  XIII 
YMAX<1U  SAVEI7    , 1 1 , ERROR ( 1 1 , PSAVE( K7 , U , 


I 


THE  COEFFICIENTS 
ORDER,  THEREFORE 


DATA    PERTST 


IN 

ONL 


/12.G 

12.0 
1  1.  ,  6.  ,  2.  ,  •  3 

DATA   MF,MTYP    /2,1 
DATA    A(2I    /    -1.0/ 
IF( ITPSC    .EO.    1     J 
!         I  TPSC    »    0 
N14=N/M0 
JRET   =2 
K2=K7/2 
K3=K7/2+l 
IF<K7/2*2.NE.K7I 
KFLAG    »    -5 
RETURN 
IRET    =    1 
KFLAG    =    1 
IF     (JSTART.LE.O) 
IF(  IQQQC.EQ.l)    GO 

BEGIN  BY  SAVING  INFO 
H  BY  THE  FACTOR  R  IF 
DEPENDENT  ON  H  MUST 
E  IS  A  COMPARISON  FO 
TO  TEST  FOR  INCREASI 
HNEW  IS  THE  STEP  SIZ 


PERTST  ARE  USED  IN  SELECTING  THE  STEP  AND 
Y  ABOUT  ONE  PERCENT  ACCURACY  IS  NEEDEO. 

»12.0,9.60,9.231,1.,1. ,1., 
,7.2,4.  102  ,1.,1.,1.,1.  , 
,1.,1.  ,1./ 
/ 

GO  TO  231 


GO  TO  1 


GO  TO  140 
TO  141 

RMATION  FOR  POSSIBLE  RESTARTS  AND  CHANGING 

THE  CALLER  HAS  CHANGED  H.   ALL  VARIABLES 
ALSO  BE  CHANGED. 

R  ERRORS  OF  THE  CURRENT  ORDER  NO.  EUP  IS 
NG  THE  ORDER,  EDWN  FOR  DECREASING  THE  ORDER. 
E  THAT  WAS  USED  ON  THE  LAST  CALL. 


100   DO  110  I  =  1,N 

HSAV( 1,1)  =  HMY(1,I » 
HSAV(2,I)  =  HMY12,II 
DO  110  J  =  1,K 

110       SAVE! J, II  =  Y(J,II 
HOLD  =  HNEW 
IF  (H.EQ.HOLDI  GO  TO  130 

120   RACUM  =  H/HOLD 


108. 


DSB19700 
DSB19800 
DSB19900 
DSB20000 
DSB20100 
DSB20200 
DSB20300 
DSB2C400 
DSB20500 
DSB20600 
DSB20700 
DSB20800 
DSB20900 
DSB21000 
DSB21100 
DSB21200 
DSB21300 
DSB21400 
DSB21500 
DSB21600 
DSB21700 
DSB218C0 
DSB21900 
DSB22000 
DSB22100 
DSB22200 
DSB22300 
DSB22400 
DSB22500 
DSB22600 
DSB22700 
DSB22800 
DSB22900 
DSB23000 
DSB23100 
DSB23200 
DSB23300 
DSB234J0 
DSB23500 
DSB23600 
DSB23700 
DSB23800 
DSB23900 
DSB24000 
DSB2<fl00 
DSB24200 
DSB24300 
DSB24400 
DSB24500 
DSB24600 
DSB24700 
DS824800 
DSB24900 
DSB25000 
DSB25100 
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IRET1    »   1 
GO   TO    750 
130      NQOLO   ■    NQ 
TOLD    «    T 
RACUM    »    1.0 

IF(    JSTART.GT.C.AND.IQQQC.NE.l)    GO   TO    250 
GO    TO    170 

140  IF    (JSTART.EQ.-ll    GO   TO    160 
C 

C      ON    THE    FIRST    CALL,    THE    ORDER    IS    SET    TO    1    AND    THE    INITIAL 

C       DERIVATIVES    ARE    CALCULATED. 

C 

NQ    »    1 

141  N3    «    N 

NMO    »    N   -    MO 

N1=N*7 

L3L    *    9*N 

N2    »   Nl    ♦    1 

N4   *    N**2 

N5    *   Nl    ♦    H 

N6   *    N5    ♦    I 

NY1    »    N*3 

NY2=NY1+1 

NY5=NY1+N 

NY6=NY5*1 

IF( JSTART.GT.O)       GO    TO    100 

JRET=1 

Ad  »    '-.5000000000000 

CALL    DIFFUN(X,Y,T, DIVER, MO, INT, N14, 1 DER , ERR ,SAVE< N2  ,1  1 , 1,  0,  0» 

EPSN   =    .167D0    *    EPS/DFLOAT(N» 

DO    144    1=1, N 

Nil    =    Nl    ♦    I 
Y(2,I  >    =    SAVE(N11,1»*H 

ERROR( II  =  EPSN 
144   CONTINUE 

CALL  DIFFUN(X,Y,T,DIVER,M0,lNT,N14,IDER,ERR,SAVE(N2,n,0,l,ll 
GO  TO  310 

142      CONTINUE 

CALL    DIVERR(PE0,N,ERR0R,ERR,HSAV(NY2,1»,HMY,PERTST(NQ,1I,K7,1,M0 
1 ,  X  I 

J  RET*  2 

DO    150    I    =*    1,N 

ERROR  (I)    »    0.  0 

NY11=NY1+I 

HMY(2,I)»HSAV(NY11,1»*H 

150   CONTINUE 
HNEW  »  H 
K  *  2 
GO  TO  100 
C 

C   REPEAT  LAST  STEP  BY  RESTORING  SAVED  INFORMATION. 
C 
160   IF  (NQ.EO.NQOLDI  JSTART  «  1 

T»TOLD 

NO  *  NQOLD 


DS825200 

DSB25300 

DSB25400 

DSB25500 

DSB25600 

DSB25700 

DSB25800 

DSB25900 

DSB26000 

DSB26100 

DSB26200 

DSB26300 

DSB26400 

DSB26500 

DSB26600 

DSB26700 

DSB26800 

DS826900 

DSB27000 

DSB27100 

DSB27200 

DSB27300 

DS827400 

DSB27500 

DSB27600 

DSB27700 

DSB27800 

DSB27900 

DSB28000 

DSB28100 

DSB28200 

DSB28300 

DSB28400 

DSB28500 

DSB28600 

DSB28700 

DSB28800 

0SB28900 

DSB29000 

DSB29100 

DSB2920C 

DSB29300 

DSB29400 

DSB29500 

DSB29600 

DSB29700 

DS829800 

DSB29900 

DSB30000 

DSB30100 

DSB30  200 

DSB30300 

DSB30400 

DSB30500 

DSB30600 
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K 
GO 


■    NQ   ♦    1 
TO    120 


SET    THE    COEFFICIENTS    THAT   DETERMINE    THE    ORDER    AND   THE    METHOD 
TYPE.       CHECK    FOR    EXCESSIVE    ORDER.       THE    LAST    TWO    STATEMENTS    OF 
THIS    SECTION      SET    IWEVAL    .GT.J     IF    PW    IS    TO    BE    RE-EVALUATED 
BECAUSE    OF    THE    ORDER    CHANGE,    AND    THEN    REPEAT    THE    INTEGRATION 
STEP    IF    IT    HAS    NOT    YET    BEEN    DONE    (IRET    *    il    OR    SKIP    TO    A    FINAL 
SCALING    BEFORE    EXIT    IF    IT    HAS    BEEN    COMPLETED    (IRET    ■    21. 


170       IF(NQ    .GT.    41    GO    TO    190 

GO   TO    (221, 222, 223, 224), NQ 

190      KFLAG    »    -2 
RETURN 


C 

C 
C 
C 

c 
c 
c 

c 
c 


THE  FOLLOWING  COEFFICIENTS  SHOULD  BE  DEFINED  TO  THE  MAXIMUM 
ACCURACY  PERMITTED  BY  THE  MACHINE.   THEY  ARE  IN  THE  ORDER  USED., 

-1/2 
-1/2,-1/2 

-1/2,-5/8,-1/8 
-1/2,-4  7/72,-1/6,-1/72 


221  A(l  )*-. 5000000000 
GO  TO  230 

222  A(l »=-. 5000000000 
A(3)  =  A(1» 

GO  TO  230 

223  AC1)=-. 5000000000 
A(3»=-.6250000C0C0 
A(4)=-. 125000000000 
GO    TO    230 

224  A(l »=-. 500C0000000 

A(3 )=-. 6527777777777778 
A(4)=-. 1666666666666667 
A( 5 »=-. 01 38888888 8888889 

230  K    =    NQ>1 
IDOUB    =    K 

231  CONTINUE 

ENQ2    =    .5/FL0AT(NQ    ♦    II 
ENQ3    =    .5/FL0AT(NQ    ♦    2» 
ENQ1    =    0.5/FL0AT(NQI 
PEPSH    =    EPS 

EUP    =     (PERTST(NQ,2J*PEPSH)**2 
E    =    (PERTST(N0,1»*PEPSHI**2 
EDWN=(PERTST(NQ,3I*PEPSH»**2 
IF    (EDWN.EQ.OI    GO    TO    780 
BND    =    EPS*ENQ3/DFL0AT(N> 
IF(  ITPSC    .EQ.    1     I    GO    TO    2 
240       IWEVAL    =    MF 

GO    TO    (    2  50    ,    680    >,IRET 


THIS    SECTION    COMPUTES    THE    PREDICTED    VALUES    BY    EFFECTIVELY 
MULTIPLYING    THE    SAVED    INFORMATION   BY    THE    PASCAL    TRIANGLE 
MATRIX. 


OSB30  700 

DSB30800 

DSB30900 

DSB31000 

DSB31100 

DSB31200 

DSB31300 

DSB31400 

DSB31500 

DSB31600 

DSB31700 

0SB31800 

DSB31900 

DSB32000 

DSB32100 

DSB32200 

0SB32300 

DSB32400 

DS832500 

DSB32600 

DSB32700 

DSB32800 

DSB32900 

DSB33000 

DS833100 

DSB33200 

DSB33300 

DSB33400 

DSB33500 

DSB33600 

DSB33700 

DSB33800 

DSB33900 

DSB34000 

DSB34100 

DSB34200 

DSB34300 

DSB34400 

DSB34500 

DSB34600 

DSB347C0 

DSB34800 

DSB34900 

DSB35000 

DSB35100 

DSB35200 

DSB35300 

DSB35400 

DSB35500 

DSB35600 

DSB35700 

DSB35800 

DSB35900 

DSB36000 

DSB36100 


c 
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c 

250   T  =  T  ♦  H  DSB36200 

DO  260  J  »  2,K  DSB36300 

DO  260  Jl  =  J,K  DSB36400 

J2  »  K  -  Jl  ♦  j  -  1  0S836500 

00  260  I  =  1,N  0SB36600 

^  ,f,j  •%n*i\^i£»tvA\»?:l\un  *  HHv,l•,,  ♦  H"Y,^•I,    lilHsH 

C  UP    TO   3   CORRECTOR    ITERATIONS    ARE   TAKEN.      CONVERGENCE    IS   TESTED  n1«l7?n? 

C  tV^Or'tIst^OnItanT."    L6"    ™N    "»    """    «    ™™™°»  "»™ 

S     ™5!!f,?.w  ?!  I^SM'Si  ^^SkV^^,*      i!?° 

?     rl    TH!^iFACTORIAL(K-ll»A(KI»,  AND  IS  THEREFORE  PROPORTIONAL       DSB37600 
C     TO  THE  ACTUAL  ERRORS  TO  THE  LOWEST  POWER  OF  H  PRESENT.  IH**K»       DSB377^S 

DO  270  I  =  1,N  DSB3780O 

270     ERRORUI  *  0. G  0SB37900 

Li=3  DSB38000 

DO  430  L  =  1.L1  ,  DSB38100 

CALL  DIFFUNU,Y,T,DIVER,M0,INT,N14,IDER,  ERR,  SAVECN2,1»,  0,0,01  ollllllo 

llrl^J,    HAS  BEEN  A  CHANGE  0F  0RDER  0R  THERE  HAS  BEEN  TROUBLE  DSB38^0 

C     WITH  CONVERGENCE,  PW  IS  RE-EVALUATED  PRIOR  TO  STARTING  THE  DSR3fl^n 

C  THF^E,Cp?RrnTE?AII0N    IN   ™E    CASE    °F    STIFF    METH3°S.       IWEvIl    IS  ollll^O 

C  THEN    SET    TO    -1    AS    AN    INDICATOR    THAT    IT    HAS    SEEN    DONE.  DSB38800 

IFdWEVAL.LT.il    GO    TO    370  nfo!oo2° 

GO    TO    310  DSB39000 

C  DSB39100 

C  SETS    THE    JACOBIAN    FOR    USE    WITH   ASYMPTOTIC    ERROR    EQUATION  DSBsJsSo 

271  00    272    1=1, K7  DSB39400 
DO    272    J=1,N  DSB39500 

272  PEQ(I,J)=PSAVE(I,J)  DSB3960O 
C  DSB39700 

C  PREPARES    TO    DECOMPOSE    THE    MATRIX    PSAVE  DSbI'wo 

IFURET.EQ.il    GO    TO    142  n?^°?00 

IFJJRET.EQ.3I    GO    TO    721  Sco^JSS 

273  R  =  A(1»*H  DSB40200 
DO  28C  I  =1,«7  DSB4D300 
K8=K2+2-I  DSB4G400 
K8=  MAX0{1,K8I  DSB40500 
K9=N*K2  +  1-I  DS6A0600 
K9=MINl'(N,K9|  DSB40700 
DO  28C  J=K8,K9  DSR4C8C0 
IF(  J  .LE.  MO  .OR.  J  .GT.  NMO.  GO  TO  275  DSB^O 
PSAVEd.JI  =  PSAVE(I,JI>R  n«2J?«n 
GO    TO    280  DSB4U00 

275      PSAVEd.JI     =    PSAVE(I,JI     *    A(ll  S?2tJ!?2 

280      CONTINUE  DSB41330 

29C    DO    300    1=1, N  DSB41400 

IF(     I    .LE.    MO       .OR.        I    .GT.     NMO  I       GO    TO    300  SIbJuOO 
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PSAVE(K3,II    »    1.0    ♦    PSAVE(K3,I) 

300      CONTINUE 

IWEVAL    »   -1 

DO    308    I    »    It  MO 

J    »    K2-I+2 

J3    =    0 

DO    303    J2    «    JtKT 

IF(DABS(PSAVE(J2, I) )    .LT.    1.0D-20)    J3    »    J3    ♦    1 
303      CONTINUE 

J    =    K7-J+1 

IFU3    .GE.    J)    PSAVE(K3,II    =    l.ODO 

J3    =    C 

Jl    =    NMO    ♦    I 

J    =    K3+M0-I 

DO    305    J2    =    ltJ 

IF(DABS(PSAVE<J2,J1 )).LT.    l.CD-20)    J3    =    J3    «•    1 
305      CONTINUE 

IF(J3    .GE.    Jl    PSAVE(K3,J1I    »    l.ODO 
308      CONTINUE 

CALL       0ECOMP<N,PSAVE,Jl,K7l 

'  CALCULATE   THE    JACOBIAN   BY    NUMERICAL    DIFFERENCING 


IF( J1.GT.0)    GO    TO    370 
WRITE<6, 11112) 
11112    FORMATC    MATRIX    PROBLEM*) 
GO    TO   440 
310  DO   320    I    =    1,N 

320  SAVE(6,I)    =   Yd, I) 

DO    340    J    =    1,N 

R    =    EPS*DMAX1(EPS,DABS(SAVE(6,J))) 
Y(l, J)    =    Yd, J)    ♦    R 
D=l./R 

CALL    DIFFUN(X,Y,T,DIVER,M0,INT,N14,IDER,ERR,SAVE<N6,l ), 0,1,1) 
DO    335    I    =    1,K7 
JI=J-K3+I 

IFUI.LE.O)  GO  TO  335 
IF(JI.GT.N)  GO  TO  340 
Il  =  K7d-I 
N12=N5+JI 
N13=N1+JI 
330   PSAVE(I1,JI)  =  (SAVE(N12,1)-SAVE(N13,1H«D 
335   CONTINUE 
340       Yd, J)  =  SAVE(6,J) 

GO  TO  271 
350   CONTINUE 
370       DO  380  I  =  1,N 
Nil  =  N5  ♦  I 
N12  =  Nl  ♦  I 
SAVE(N11,1)   =Y(2,I)   -   SAVE(N12,1)*H 
IF(I  .LE.  MO   .OR.   I  .GT.  NMO)  SAVE<Nll,l)  =  -SAVE(N12,1) 
380   CONTINUE 
C 

C     COMPLETE  THE  GAUSSIAN  ELIMINATION 
C 


OSB41700 

DSB41800 

DSB41900 

DSB42000 

DSB42100 

DSB42200 

DSB42300 

DSB42400 

DSB42500 

DSB42600 

DSB42  700 

DSB42800 

DSB42900; 

DSB43C00 

DSB43100j 

DSB43200! 

DSB43300 

DSB43400 

DSB43500 

DSB43600 

DSB43700 

DSB43800 

DSB43900 

DSB44000 

DSB44100 

DSB44200 

DSB44300 

DSB44400 

DSB44500 

DSB44600 

DSB44700 

DSB44800 

DSB44900 

DSB45000 

DSB45100 

DSB45200 

DSB45300 

DSB45400, 

DSB45500 

DSB45600 

DSB45700 

DSB45800 

DSB45900 

DSB46000 

DSB46100 

DSB46200 

DSB46300 

DSB46400 

DSB46500 

DSB46600 

DSB46700 

DSB46800 

DSB46900 

DSB47000 

DSB47100 
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CALL    S0LVE(N,PSAVE,SAVE<N5-H,ll,SAVE,K7,7,6l 
410  NT    «    N 

;  PERFORM   THE    CORRECTION    STEP 


ltN 

Yd, II    ♦ 
Y(2,I)    - 


A(1>*SAVE(6,I I 
SAVE(6,II 


42  0 


NT 


NT   -    1 


42  5 


42  7 


430 


00  420  I 
Y(1,I) 
Y(2,I> 

ERRORUI    »    ERRORm    ♦    SAVE(6,I) 
IF       (DABS(SAVE<6,II  >.LE.  ( BND*YMAX < 1 1  I  I 
CONTINUE 

CALL  DIVERR(PEQ,N,ERROR,ERR,HSAV(NY2,1I,HMY,PERTST(NQ,1I,K7,0,MO 
lfXI 
DO  4251  =  1, N 
NY11=NY5+I 
NY12=NY1+I 

HSAV  (NYll,l»=HMY(2,n-HSAV(NY12,ll*H 

IF(  I    ,LE.    MO      .OR.       I    .GT.    NMO)    HSAV(NY11,1|    *    -    HSAV(NY12,1» 
CONTINUE 

SOLVE    THE    ASYMPTOTIC    ERROR    EQUATION 

CALL    SOLVE(N,PSAVE,HSAV(NY5*l,ll,HSAV,K7,3,3) 
DO    427    I»1,N 

HMYt 1 , I »  =  HMY ( 1 , II +A ( 1 l*HSAV( 3, I  I 
HMY(2,I)=HMY(2,II-  HSAV<3,II 

IF    (NT.LE.OI    GO    TO   490 

CONTINUE 


THE    CORRECTOR    ITERATION    FAILEO    TO    CONVERGE    IN    3 

POSSIBILITIES    ARE    CHECKED    FOR,        IF    H    IS    ALREADY 

THIS    IS    EITHER    ADAMS    METHOD    OR    THE    STIFF    METHOD 

MATRIX    PW    HAS    ALREADY    BEEN    RE-EVALUATED,    A    NO    CONVERGENCE    EXIT 

IS    TAKEN.    OTHERWISE    THE    MATRIX    PW    IS    RE-EVALUATED    AND/OR    THE 

STEP    IS    REDUCED    TO    TRY    AND    GET    CONVERGENCE, 


TRIES,       VARIOUS 

HMIN    AND 

IN    WHICH    THE 


440   T=TOLD 

IF     UH.LE.  (HMIN*1.  00001  ) ) . AND. ( ( I WEVAL    -    MTYP  I.  LT.-l  I  I 
IF     HMF.EQ.OI.OR.  (  I  WEVAL.  NE  ,0  I  »    RACUM    =    RACUM*0.25DO 
IWEVAL    =    MF 
IRET1    «    2 
GO    TO    750 
KFLAG    *    -3 
DO    480    I    » 
HMY(lfl)     * 
HMYC2.II    = 
DO    480    J 
Y(J,II 
H    =    HOLD 
NQ    =    NQOLD 
J  ST  ART    «    NQ 
RETURN 


GO    TO    460 


460 
470 


480 


i,N 

HSAV(1,II 
HSAV(2,I I 
=    i,K 
=    SAVEIJtll 


THE    CORRECTOR    CONVERGED    ANO    CONTROL     IS 
IF    THE    ERROR    TEST    IS    O.K.,       AND    TO    540 


PASSED    TO    STATEMENT 
OTHERWISE, 


520 


DSB47200 
DSB47300 
DSB47400 
DSB47500 
DSB47600 
DSB47700 
DSB47800 
DSB47900 
DSB48000 
DSB48100 
DSB48200 
DSB48300 
DSB48400 
DSB48500 
DSB48600 
DSB48700 
DSB48800 
DSB48900 
DSB49000 
DSB49100 
DSB49200 
DSB49300 
DSB49400 
DSB49500 
DSB49600 
DSB49700 
DSB49800 
DSB49900 
DSB50000 
DSB50130 
DSB502OO 
DSB50300 
DSB50400 
DSB50500 
. DSB50600 
DSB50700 
DSB508D0 
DSB50900 
DSB51000 
DSB51100 
DSB51200 
DSB51300 
DSB51400 
DSB515C0 
DSB51600 
DSB51700 
DSB51800 
DS851900 
DSB52000 
DSB52100 
DSB52200 
DSB52300 
DSB52400 
DSB52500 
DSB52600 
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C   IF  THE  STEP  IS  O.K.  IT  IS  ACCEPTED.  IF  I DOUB  HAS  BEEN  REDUCED 
C   TO  ONE,   A  TEST  IS  MADE  TO  SEE  IF  THE  STEP  CAN  BE  INCREASED 
C   AT  THE  CURRENT  ORDER  OR  BY  GOING  TO  ONE  HIGHER  OR  ONE  LOWER. 
C   SUCH  A  CHANGE  IS  ONLY  MADE  IF  THE  STEP  CAN  BE  INCREASED  BY  AT 
C   LEAST  1.1.   IF  NO  CHANGE  IS  POSSIBLE  I DOUB  IS  SET  TO  10  TO 
C   PREVENT  FURTHER  TESTING  FOR  10  STEPS. 

C   IF  A  CHANGE  IS  POSSIBLE,  IT  IS  MADE  AND  IDOUB  IS  SET  TO 
C   NQ  ♦  1     TO  PREVENT  FURTHER  TESING  FOR  THAT  NUMBER  OF  STEPS. 
C   IF  THE  ERROR  WAS  TOO  LARGE,  THE  OPTIMUM  STEP  SIZE  FOR  THIS  OR 
C   LOWER  ORDER  IS  COMPUTED,  AND  THE  STEP  RETRIED.   IF  IT  SHOULD 
C   FAIL  TWICE  MORE  IT  IS  AN  INDICATION  THAT  THE  DERIVATIVES  THAT 
C   HAVE  ACCUMULATED  IN  THE  Y  ARRAY  HAVE  ERRORS  OF  THE  WRONG  ORDER 
C   SO  THE  FIRST  DERIVATIVES  ARE  RECOMPUTED  AND  THE  ORDER  IS  SET 
C   TO  1. 
C 
490   D  =  0.0 

DO  5C0  I  »  ltN 
500     D  =  D  ♦  (ERR0R(I)/YMAX( 1)1**2 

IWEVAL  *  0 

IF  (D.GT.EI  GO  TO  540 

IF  (K.LT.3)  GO  TO  520 
C 
C     STEP  SUCCEEDED,   UPDATE  REMAINING  VECTORS 


510 
520 


3,K 

*    1,N 

«    Y(J,I)    ♦    A(J)*ERROR(I) 


DO    510    J    - 
DO   510    I 
Y( J, I  I 
KFLAG    *    +1 
HNEW   =    H 

IF    ( IDOUB. LE. II    GO    TO   550 
IDOUB    =    IDOUB    -    1 
IF    (IDOUB.GT. 1)    GO    TO   700 
DO    53C    I    =    1,N 
530  SAVEI7    ,1)    =»    ERROR(I) 

GO    TO    700 

ERROR    TOO   LARGE,    REDUCE    FAILURE    FLAG    AND    TRY    AGAIN 

540      KFLAG    =    KFLAG   -    2 

IF    (H.LE.  (HMINM.OOOCl)  )    GO    TO   740 

T    =    TOLD 

IF    (KFLAG. LE. -11)       GO   TO    720 

CHECKS    TO    SEE    IF    CHANGES    SHOULD    BE    MADE    IN    STEP    SIZE    OR    ORDER 

550      PR2    =    (D/E)**ENQ2*1.2 
PR3    =    l.E+20 

IF    UNQ.GE.MAXDER). OR. (KFLAG. LE.-l)  )    GO    TO    570 

D    =    O.C 

DO  560  I  =  1,N 
560     D  =  D  ♦  ((ERROR(I)  -  SAVE(7  , I )  )/YMAX( I ) 1**2 

PR3  =  (D/EUP)**ENQ3*1.4 
570   PR1  =  l.E+20 

IF  (NQ.LE.l)  GO  TO  590 

D  =  0.0 


DSB52700 

DSB52800 

DSB52900 

DSB53000 

DSB53100 

DSB53200 

DSB53300 

DSB53400 

DSB53500 

DSB53600 

DSB53700 

DSB53800 

DSB53900 

DSB54000 

DSB54100 

DSB54200 

DSB54300 

DSB54430 

DSB 54500 

DSB54600 

DSB 54700 

DSB54800 

DSB54900 

DSB55000 

DSB55100 

DSB55200 

DSB55300 

DSB55400 

DSB55500 

DSB55600 

DSB55700 

DSB55800 

DSB55900 

DSB56000 

DSB56100 

DSB56200 

DSB56300 

DSB56400 

DSB56500 

DSB56600 

DSB56700 

DSB56830 

DSB56900 

DSB57000 

DSB57100 

DSB57230 

DSB57300 

DSB57400 

DSB57500 

DSB57600 

DSB57700 

DSB57800 

DSB57900 

DSB58000 

DSB58100 
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00  580  I  *  i,N 
580     0  -  0  ♦  (Y(K,I l/YMAX(I>l**2 

PR1  *  (D/EDWN>»*ENQ1*1.3 
590   CONTINUE 

IF  (PR2.LE.PR3)  GO  TO  650 

IF  (PR3.1T.PRU  GO  TO  660 
600   R  »  1.0/AMAXKPR1  ,l.E-4> 

NEWQ  «  NQ  -  1 
610   ID0U8  »  10 

IF  ((KFLAG.EQ.ll. AND.(R.LT.  (l.lil )  GO  TO  700 

IF  (NEWQ.LE.NQ)  GO  TO  630 


COMPUTE  AN  ADDITIONAL  SCALED  DERIVATIVE  IF  ORDER  IS  INCREASED 

DO  620  I  =  1,N 

620     Y(NEWQ*1,II  «  ERRORd  »*A(K)  /DFLOAT(  Kl 
630   K  =  NEWQ  ♦  1 

IF  (KFLAG.EQ.ll  GO  TO  670 
RACUM  =  RACUM*R 
IRET1  »  3 
GO  TO  750 
640   IF  (NEWQ.  EQ.NQ)  GO  TO  250 
NQ  *  NEWQ 
GO  TO  170 
650   IF  (PR2.GT.PR1)  GO  TO  600 
NEWQ  *  NQ 

R  =  1.0/AMAXKPR2,l.E-4) 
GO  TO  610 
660   R  =*  1.0/AMAX1(PR3,1.E-4I 
NEWQ  »  NQ  ♦  1 
GO  TO  610 
670   IRET  »  2 

R  =  DMIN1(R,HMAX/DABS(H)| 
H  »  H^R 
HNEW  =  H 

IF  (NQ.  EQ.NEWQ)  GO  TO  680 
NQ  =  NEWQ 
GO  TO  170 
680   Ri  =  1.0 

DO  690  J  =  2,K 
Rl  =  R1*R 
DO  690  I  =  1,N 
IF(J  .EQ.  2)   HMY(J,IJ  »  HMY(J,I»  *R1 
690       Y(J,I»  =  Y(J,n»Rl 

IDOUB  =  K 
700   DO  710  I  *  1,N 

YMAX(I)  =  DMAX1(YMAX(II,DABS(Y(1,I))J 
710   ER0RA(n=(HMY(l,  I  l-HSAVI  1  ,11) /YMAXdl 
JSTART  »  NQ 
RETURN 
720   IF  (NQ.EQ.1J  GO  TO  780 

CALL  DIFFUN(X,Y,T, DIVER,  MO,  I  NT, N14, 1 DER, ERR , SAVE(  N2, 1 » , 1,  0 , 0 1 
J  K  c  1*3 

EPSN  =  .16700  ♦  EPS/OFLOAT(NI 
DO  722  1=1, N 


0S858200 

0SB58300 

DSB58400 

DSB58500 

0SB58600 

DSB58700 

DSB58800 

0SB58900 

DSB59000 

0SB59100 

DSB59200 

DSB59300 

DSB59400 

DSB59500 

DSB59600 

DSB59700 

DSB59800 

DSB59900 

DSB60000 

DSB60100 

OSB60200 

DSB60300 

DSB63400 

DSB60500 

0SB60600 

DSB60700 

DSB60800 

DSB6C900 

DSB61000 

DSB61100 

0SB61200 

DSB61300 

DSB61400 

DSB61500 

DSB61600 

0SB61700 

DSB61800 

DSB61900 

DSB62000 

DSB62100 

DSB62200 

DSB62300 

DSB62400 

DSB62500 

DSB62600 

DSB62700 

DSB628G0 

DSB62900 

DSB63000 

DSB63100 

DSB63200 

DSB63300 

DSB63400 

DSB63500 

DSB63600 
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722 
721 


730 


740 


C 

C       TH 
C      TO 
C 
750 


760 


770 


780 


R  .«    H/HOLO 

Yd, II    «    SAVEd.I  I 

Nil    *    Nl    ♦    I 

SAVE(2,I)    ■    HOLD*SAVE(NU,l» 
Y(2,I  I    =    SAVE(2,I  I    *    R 

ERRORd  I    «    EPSN 
CONTINUE 

CALL    DIFFUN(X,Y,T,DIVER,MQ, INT,N14,IDER,ERR,SAVE(N2,ll,0,i,il 

GO    TO    310 

CONTINUE 

JRET-2 

CALL    DIVERR(PEQ,N,ERR0R,ERR,HSAV(NY2,1),HMY,PERTST(NQ,1I,K7,1,M0 
1,XI 
R    =    H/HOLD 
DO    73C    I    *    1,N 
ERRORd  I    »    0,0 
HMY(l,II=HSAVd,I  I 
NY11=NY1*I 

HSAV(2,I)*H0LD*HSAV(NY11,1> 
HMY(2,I l=HSAV(2,I l*R 
CONTINUE 
NQ  *  1 
KFLAG  »  1 
GO  TO  170 
KFLAG  *  -1 
HNEW  *  H 
JSTART  *  NQ 
RETURN 

IS  SECTION  SCALES  ALL  VARIABLES  CONNECTED  WITH  H  AND  RETURNS 
THE  ENTERING  SECTION. 

RACUM  =  DMAX1 ( DABS ( HMIN/HOLD I , RACUM ) 
RACUM  =  DMIN1  (RACUM, DABS( HMAX/HOLDII 
Rl  =  1.0 
DO  760  J  =  2,K 

Rl  =  R1*RACUM 

DO  760  I  =  1,N 
IF(  J  .EG.  21  HMY(J,II  =  HSAV(J,II  *  Rl 

Y(J,  I  I  =  SAVE(J,I I*R1 
H  =  HOLD*RACUM 
DO  770  I  =  1,N 

HMYd,II=HSAVd,I  I 

Yd, II  =  SAVEd  ,1  I 


IDOUB 
GO  TO 
KFLAG 
GO  TO 
END 


=  K 
(  130 
=  -4 

470 


,  2  50  ,  640  »,IRET1 


DSB63700 

DSB63800 

DSB63900 

DSB  640001 

DSB64100 

DSB64200 

DSB6430C 

DSB64400 

DSB64500 

DSB 64600 

DSB64700 

DSB64800 

DSB64900! 

DSB65000 

DSB651001 

DSB65200 

DSB65300 

DSB65400 

DSB65500 

DSB65600 

DSB65700 

DSB65800 

DSB65900 

DSB66000 

DSB66100 

DSB66200 

DSB66300! 

DSB66400J 

DSB66500  i 

DSB66600 

DSB66700 

DSB66800 

DSB66900 

DSB67000 

DSB67100 

DSB67200': 

DSB67300" 

DSB67400'1 

DSB67500 

DSB67600 

DSB67700 

DSB67800 

DSB67900 

DSB68000 

DSB68100 

DSB68200 

DS868300 

DSB68400 

DSB 68500 


SUBROUTINE  DI VERR (PEQ,N, ERROR  ,ERR,DY, Y, PERTST ,K7,M,M0 ,X) 


DVR10000 
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c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


THE  PURPOSE  OF  THIS  SUBPROGRAM  IS  TO  CALCULATE  THE  TIME 
DERIVATIVE  OF  THE  ASSYMPTOTIC  ERROR  EQUATION 

THE  VARIABLES  IN  THIS  SUBPROGRAM  ARE 


PEQ 
DY 

N 

ERROR 

ERR 

Y 

PERT  ST 

K7 


PARTIALS  OF  THE   DIFFERENTIAL  EQUATIONS 

WRT  THE  VARIABLES 

THE  TIME  DERIVATIVE  VECTOR  OF  THE  ASYMPTOTIC 

ERROR  EQUATION 

NUMBER  OF  EQUATIONS 

THE  TIME  BASED  TRUNCATION  ERROR  VECTOR 

THE  SPACE  BASED  TRUNCATION  ERROR  VECTOR 

VALUE  VECTOR  OF  THE  ASYMPTOTIC  ERROR 

CONSTANTS  ASSOCIATED  WITH  THE  ORDER 

THE  NUMBER  OF  DIAGONALS  IN  THE  JACOBIAN 


IMPLICIT  REAL  *  8  IA-H.O-ZI 

REAL  *  4  PEPTST 

COMMON  /DVR/  DVYC4) ,COEF ( 4, 3 ,3 ) , ALPHA ( 2,4, 3 > , DI V( 16, 5 1 

DIMENSION  DY(1» ,P EQ < K7, 1 ) ,E RROR(l) ,ERR<1>  ,Y(2,1I 


DIMENSION  XT(5I ,XS< 5 1 ,YT< 3, 5 1 ,VSI  3,  5  1 

DIMENSION  D(4),ER<4> 

DIMENSION  X(l» 

K8  *  K7/2  +1 

K2  »  K8  ♦  I 

K3  *  K8  ♦  N 

DO  10  IA  »  1,N 

Ni  =  MAXO(K2-IA,l  ) 

N2  =  MIN0U3-IA.K7I 

Kl  =  IA  -  K8 

c 

c 
c 

DY(IA)»-     (ERR(IA)I   ♦   ( ERROR! I A J/PERTSTI 

MULTIPLIES  THE  APPROPRIATE  VARIABLES  BY  THE  PARTIA 

DO   5  IB  =  Nl  ,  N2 

DY(IA)  *  PEQ(IB.IA»*YU,K1*IB>  ♦  DY(IA» 

5 

CONTINUE 

10 

c 
c 
c 

CONTINUE 

THIS  SECTION  DETERMINES  THE  VALUES  AT  THE  BOUNDARY 

IF(M  .EQ.  0)  RETURN 

L=N/MO 

NO  *  i 

xi  =  xm 

XN  »  X(LI 

DO  115   J  *  1,M0 

DO   110  I  »  1,5 

11  »  I  ♦  1 

12  *  L  -  I 

XT(  I)  =  Xdll 

XS(  II  ■  X  (  I  2  1 

11  =  I*MO*J 

YT(1,II  *  DY(  H  1 

12  *  M0*(I2-1)  ♦  J 

DVR10100 

DVR1D200 

DVR10300 

DVR10400 

DVR10500 

DVR10600 

DVR13700 

DVR10800 

DVR10900 

DVR11000 

DVR11100 

DVR11200 

DVR11300 

DVR11400 

DVR11500 

DVRH600 

DVR11700 

DVR11800 

DVR11900 

DVR12000 

DVR12100 

DVR12200 

DVR12300 

DVR12400 

DVR12500 

DVR1260O 

DVR12700 

DVR128C0 

DVR12900 

DVR13000 

DVR13100 

DVR13200 

DVR1330O 

DVR13400 

DVR1350O 

DVR13600 

DVR13700 

DVR13800 

DVR13900 

DVR14000 

DVR14100 

DVR14200 

DVR14300 

DVR144Q0 

DVR14500 

DVR14600 

DVR147C0 

DVR14800 

DVR14900 

DVR15000 

DVR15100 

DVR15200 

DVR15300 

DVR15400 

DVR15500 
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YSCl.II    -   DYU2I  .UB1C1AA 

110      CONTINUE  DVR15600 

CALL    DERIVCXT,YT,IDER,NO,0,l,C,Na,Xl,D,XII  Dwlslnn 

CALL    DERIV(XS,YS,IDER,NO,OU.O,NO.XN,ER,XII  dvrisq™ 

dy(jj  »  Ddi  n»„:;;00 

n  *  mo*(l-ii*j  RXSJ!?00 

DYCIll  »  ERC1I  2X?}6100 

115   CONTINUE  DVR16200 

RETURN  DVR16300 

EN0  DVR16400 

0VR16500 


119 


APPENDIX  B 
EXAMPLE  1 
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903 


30 


50 


IMPLICIT  REAL  *81A-H,Q-Z» 

DIMENSION  Y(2,513),YTST<513»,M{8),E<513»,F{513I,TTST(8I 

DATA  M/3, 7, 15, 31, 63, 127, 255, 511/ 

DATA  TTST/.125,.25,.5,1.,2.  ,4,  ,8.  ,16.  / 

RP I =3. 141 592653  589793 

QPI»RPI+RPl 


IP1 
IP2 
IP3 
IP4 
IP5 
IP6 
IP7 


100 


MAIN    LOOP    FOR    DIFFERENT    STEP    SIZES 

DO    500    IA    «    2,7 

RMAX    *    0,0 

ID    =   1 

KA»1 

KB»2 

T=O.D 

NA=M( IA I 

WRITE(6,903I    NA 

F0RMATP1SI30I 

NB=NA+1 

NUMM=«0 

NUMD=0 

NC=NA+2 

HX=1./DFL0AT(NB» 

HT=HX 

ZPI    =    HX    *    RPI 

INITILIZE 

DO    30     IB    =    1,NC 

Y(KA, IB>=DSIN(ZPI*DFLOAT< IB-ll) 

A=.5*HT/(HX*HX> 

B=1.+2.*A 

V=1.-2.*A 

C=A 

NUMM=NUMM+4 

NUMD    =    NUMD    ♦    1 

E(ll=C.O 

F(1I  =  0.0 

ONE  STEP  OF  THE  INTEGRATION  -  SOLVE  SYSTEM  OF  EQUATIONS 

Y(KA,1)  *  0.0 

DO  100  IB  =  2,NB 

CPT=B-C*E( IB-ll 

E(IB»=         A/CPT 

D=  A*Y(KA,IB-1)+V*Y(KA,IB)*A*Y(KA,IB*1I 

F(IB)  =  (D+C*F(  IB-1  H/CPT 


EXA10000 

EXAniOC 

EXA10200 

EXA10300 

EXA10400 

EXA10500 

EXA10600 

EXA10700 

EXA13800 

EXA10900 

EXA11000 

EXAiliOO 

EXA11200 

EXA11300 

EXA11400 

EXA11500 

EXA11600 

EXA11700 

EXA11800 

EXA11900 

EXA120CO 

EXA12100 

EXA12200 

EXA12300 

EXA12400 

EXA12500 

EXA12600 

EXA12700 

EXA12800 

EXA12900 

EXA13000 

EXA13100 

EXA13200 

EXA13300 

EXA13400 

EXA13500 

EXA13600 

EXA13700 

EXA13800 

EXA13900 

EXA14000 

EXA14100 

EXA14200 

EXA14300 

EXAU400 

EXA14500 

EXA14600 

EXA14700 

EXA14800 

EXA14900 

EXA15000 

EXA15100 

EXA15200 

EXA15300 
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Y(KB,NC»»0.0 
00    2DC    IB«1,NA 
IC*NC-IB 
200      Y(KB,ICI»E<  IC»»Y(KB,  IC  +  ll+FUCI 
NUMM»NUMM*6*NA 
NUMD»NUMD*2*NA 
KA=3-KA 
KB*  3-KB 
T=T*HT 
IF(    T      .LT.       TTSTCIOII    GO  TO    50 

OUTPUT  SECTION 

00  300  IB=1,NA 

VAL=DEXP(-QPI*T»*DSIN<ZPI*DFLOAT(IBM 

YTST(IB)  =  Y(KA,IB  +  1)  -  VAL 

IF(RMAX.LT.  OABS(YTST(IBM»   RMAX  «  DABS(  YTST  ( I  Bll 

CONTINUE 

WRITE  (6,9C0»  T,HX,A,B 

WRITE  (6,904)  Y(KA,IPl*l),Y(KA,IP2+ll,Y<KA,IP3+l)fY(KAfIP4+l), 
1  Y(KA,IP5  +  l),Y(KA,IP6*l)«Y(KAvIP7  +  ll 

WRITE  (6,  904  I    YTST(IP1),YTST(IP2),YTST(IP3),YTSTCIP4),YTST<IP5). 
1YTST( IP6),YTST< IP7)  »*r3i, 

F0RMAT(1X,5D26.15) 

WRITE<6,901)    NUMM,NUMD 

F0RMAT(1X,3I30) 

WRITE(6,902) 

FORMAT!*     •) 

FORMATdX, 7018.91 

ID    =    ID    ♦    1 

GO    TO    50 

.001)         CALL    EXIT 


300 


900 
901 

902 
904 


500 


IF(T    .LT. 

IF(    RMAX    , 
IP1 


I  PI 
IP2 
IP3 
IP4 
IP5 
IP6 
IP7 


16) 
LT. 
*  2 


IP2 
IP3 
IP4 
IP5 
IP6 
IP7 


CONTINUE 
CALL  EXIT 

END 


EXA15400 

EXA15500 

EXA15600 

EXA15700 

EXA15800 

EXA15900 

EXA16000 

EXA16100 

EXA16200 

EXA16300 

EXA16400 

EXA16500 

EXA16600 

EXA16700 

EXA16800 

EXA16900 

EXA17000 

EXA17100 

EXA17200 

EXA17300 

EXA17400 

EXA17500 

EXA17600 

EXA17700 

EXA17800 

EXA17900 

EXA18000 

EXA18100 

EXA18200 

EXA18300 

EXA18400 

EXA18500 

EXA18600 

EXA18700 

EXA18800 

EXA18900 

EXA19000 

EXA19100 

EXA19200 

EXA19300 

EXA19400 

EXA19500 
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APPENDIX  C 
EXAMPLE  2 


123 


903 


30 


IMPLICIT  REAL  *8<A-H,Q-ZI 

DIMENSION  YC2,5l3»,YTST(5l3l,M(8l,E(513l,F(5l31,TTSTIfll 
DATA  M/3, 7, 15, 31, 63, 127, 255, 511/  '   *  * 8* 

DATA  TTST/.125,.25,.5,l.,2.  ,4.,8.  ,16./ 
RP 1*3, 1*1592653 589793 
QPI»RPI*RPI 


IP1 
IP2 
IP3 
IP4 
IP5 
IP6 
IP7 


1 
2 
3 
4 
5 
6 
7 


50 


MAIN  LOOP  FOR  DIFFERENT  STEP  SIZES 

DO  50C  IA  ■  2,7 

RMAX  =0.0 

ID  »  1 

KA»1 

KB»2 

T=0.0 

NA=M( IA) 

WRITE(6,903»  NA 

F0RMAT(*1'VI30( 

NB=NA+1 

NUMM=0 

NUMD=*0 

NC=NA+2 

HX=1./DFL0AT(NB> 

HT  =  HX 

ZPI  =  HX  *  RPI 

INITILIZE 

DO  33  IB  «=  1,NC 

D1K  =  'hT/HX  OEXP("rtX*OFLOAT{IB-1,,*DSIN«^pI*DFLOAT(IB-in 

A  =  . 5D0*D1/HX 

B  =  1.00  ♦  2.D0  *A 

C  =  A 
Al  *  A  ♦  Dl 
A2  =  l.D0-2.D0*A+HT 
A3  =  A-Dl 
NUMM=NUMM+4 
NUMD  =  NUMD  ♦  1 
E(1»  =  0.0 
F(l »=0.0 

ONE  STEP  OF  THE  INTEGRATION  -  SOLVE  SYSTEM  OF  EQUATIONS 

Y(KA,1>  =  0.0 
DO  IOC  IB  *  2,NB 
CPT*B-C*E( IB-ll 


EXB10000 

EXB101DO 

EXB10200 

EXB10300 

EXB10400 

EXB10500 

EXB106OO 

EXB10700 

EXB10800 

EXB1C9G0 

EXB11000 

EXB11100 

EX811200 

EXB11300 

EXB11400 

EXB11500 

EXB11630 

EXB11700 

EXB11800 

EXB11900 

EXB12000 

EXB12100 

EXB12200 

EXB12300 

EXB12400 

EXB12500 

EXB12600 

EXB12700 

EXB12800 

EXB12900 

EXB13000 

EXB13100 

EXB13200 

EXB13300 

EXB13400 

EX813500 

EXB13600 

EXB13700 

EXB13800 

EXB13900 

EX814000 

EXB14100 

EXB14200 

EXB14300 

EXB14400 

EXB14500 

EXB14600 

EXB147C0 

EXB14800 

EXB14900 

EXB15000 

EXB15100 

EXB15200 

EXB15300 
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130 


200 


♦    A1*Y(KA,IB*1I 


900 

901 

902 
904 


500 


E  f I B)«  A/CPT 

D    »    A3    *    Y(KA,IB-1)    ♦    A2*Y(KA,IB» 
F( IB»*(D+C*F(IB-ll»/CPT 

Y(KB,NC)*0.0 

DO    230    IB«1,NA 

IC=NC-IB 

Y(KBfIC»»E{  IC»*Y(KB,  IC+1KFUCI 

NUMM=NUMM*6*NA 

NUMD=NUMD*2*NA 

KA=3-KA 

KB=3-KB 

T=T+HT 

IF(    T      .LT.       TTST(ID) I    GO    TO    50 

OUTPUT    SECTION 


WRITE    (6»900»    T,HX,A,8 

WRITE (6, 904>    Y(KA,IP1+1»,Y(KA,IP2*1I,Y(KA,IP3+1»,Y(KA, IP4+1J, 
1    Y(KA,IP5  +  1» ,Y(KA,IP6+1>,Y(KA,IP7+1I 

F0RMAT(1X,5D26.15) 

WRITE(6,901)     NUMM,NUMD 

F0RMAT(1X,3I30> 

WRITE(6,902» 

FORMAT!*     •» 

F0RMAT(1X,7D18.9I 

ID    =    ID    ♦    1 

GO    TO    50 


IF(T 
I  PI  = 
IP2  = 
IP3  = 
IP4  = 
IP5  = 
IP6  '- 
IP7    = 


LT, 
IP1 
IP2 
IP3 
IP4 
IP5 
IP6 
IP7 


161 
*  2 


CONTINUE 
CALL  EXIT 
END 


EXB15400 

EXB15500 

EXB15600 

EXB15700 

EXB15800 

EXB15900 

EXB16000 

EXB16100 

EXB16200 

EXB16300 

EXB16400 

EXB16500 

EXB16600 

EXB16700 

EXB16800 

EXB16900 

EXB17000 

EXB17100 

EXB17200 

EXB17300 

EXB17400 

EXB17500 

EXB17600 

EXB17700 

EXB17800 

EXB17900 

EXB18000 

EX818100 

EXB18200 

EXB18300 

EXB18400 

EXB18500 

EXB18600 

EXB18700 

EXB18800 

EXB18900 

EXB19000 
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APPENDIX  D 
EXAMPLE  3 


126 


IMPLICIT    REAL   *8(A-H,0-ZI  EXC10000 

DIMENSION    TTSTI8I  EXCIOIOO 

DIMENSION    Y(2,102  8>,M(11),YTST(1028» , RMATI 3, 1028 i ,B { 10281  EXC 10200 

DATA    M/3, 7, 15, 31, 63, 127, 256, 511, 1023, 2045, 4095/  EXC10300 

DATA    TTST/. 125, .25, .5,1. ,2.  .4,, 8., 16./  EXC10400 

CALL    UNOERZI»OFF»l  EXC10500 

RK1    »    1.    DO  EXC10600 

RK2    »    l.DO  EXC10700 

RPI*3.  141592653589793  EXC10800 

IP1    *    1  EXC10900 

IP2    *    2  EXC11000 

IP3    »    3  EXC11100 

IP4    =    4  EXC11200 

Ip5    *    5  EXC11300 

IP6    *    6  EXC1140C 

IP7    =    7  EXC11500 

C  EXC11600 

C            MAIN    LOOP    FOR    DIFFERENT    STEP    SIZES  EXC11700 

C  EXC11800 

DO    530    IA    *    2,7  EXC11900 

RMAX    «   0.0  EXC12000 

ID    *    1  EXC12100 

K**l  EXC12200 

*B»2  EXC12300 

T*0.0  EXC12400 

NA=M(IA)  EXC12500 

WRITE(6,903>    NA  EXC12600 

903       F0RMAT(«1«,I30I  EXC12700 

NB=NA*1  EXC12800 

NUMM=0  EXC12900 

NUMD=0  EXC13000 

NC=NA+2  EXC13100 

HX=1./DFL0AT(NBI  EXC13200 

QPI    =    RPI*HX  EXC13300 

HT=HX  EXC13400 

c  EXC13500 

C            INITILIZE  EXC13600 

c  EXC13700 

DO    3D    IB   =    i,NC  EXC13800 

V    =    DSIN(QPI*DFL0AT(IB-1))  EXC13900 

Y(KA,IB»    =    V*V  EXC14000 

EXC14100 

SET    UP    VALUES    FOR    MATRIX  EXC14200 

EXC 14300 

R    =    HT/(HX*HX)  EXC14400 

Bl    =    l.ODO    ♦    HX  EXC14500 

B2    =    l.ODO    ♦    HX  EXC14600 

TR    =    2. DO    *    R  EXC14700 

TRB1    =    TR*B1  EXC14800 

TRB2    =    TR*B2  EXC14900 

RNATflfll    »    -R  EXC15000 

RMAT(2,1»    =    2.00    ♦    TRB1  EXC15100 

RMAT<3,1»    *    -TR  EXC15200 

RMAT(1,NC»    »    -TR  EXC15300 


127 


RMAT(2,NC»  »  2. DO  ♦  TRB2 
RMATOtNCI-  -R 
RMAIN  *  2.D0  ♦  TR 
RMANE  »  2.00  -TR 


10 


DO  ID  I  « 
RMAT(1,II 
RMAT(2,I) 
RMAT(3,I) 
CONTINUE 


2,N8 
»    -R 
»    RHAIN 
«   -R 


DECOMPOSE    THE    MATRIX 

CALL    DEC0MP<NC,RMAT,J5,3I 

NUMM    =    5+NA 

NUMO   =    2  +  NA 

IF    <    J5    .LE.    0    1      GO    TO    200 


C  ONE    STEP    OF    THE    INTEGRATION    -    SOLVE    SYSTEM    OF    EQUATIONS 

50       B(l»    »    (2. DO    -    TRB1I    *Y(KA,1)    ♦    TR    *Y(KA,2> 
B(NC)    *    TR*YIKA,NBI     ♦    ( 2. D0-TRB2 I *Y( KA, NC  ) 
C 

C  SET    UP    THE    RIGHT    SIDE    OF    EQUATION 

C 

DO    23    I    =    2,NB 

B(I)    =    R*Y(KA,I-1I    ♦    RMANE*Y(KA,I I    ♦    R*Y<KA,I+1) 
20      CONTINUE 

CALL    S0LVE(NC,RMAT,B,Y,3, 2,KBI 
NUMM    =    NUMM    ♦    2*NC 
NUMD   =    NUMO    ♦    NC 
KA=3-KA 
KB=3-KB 
T=T+HT 

IF(    T      .LT.       TTSTUDII    GO    TO    53 
Z 
C  OUTPUT    SECTION 

c 

DO    3:0    IB=1,NA 

YTST(IB)     =    Y(KA,IB+1I    -    Y(KB,IB*1) 

IFIRMAX.LT.  DABSiYTSTUB)))   RMAX  *  DABS(  YTST  ( I  B>) 
3C0   CONTINUE 

WRITE  (6,9001  T,HX 

WRITE (6 ,904)  Y( KA , I Pl  +  1 ) , Y( KA, I P2*l ) , Y( KA , IP3  +  1 ) , Y( KA , I P4  +  1 I . 
1  Y(KA,IP54l),Y(KA,IP6*l),Y(KA,IP7U) 

i??STU(p6nY!sT!lJ7!Pn,YTST(IP2,,YTST(IP3,,YTST(IPM,YTST(IP5,» 

900  F0RMAT(1X,5D26. 151 
WRITE(6,90l)  NUMM,NUMO 

901  F0RMAT(1X,3I30) 
WRI TE(6,902» 

902  FORMAT! •  •) 

9C4   F0RMAT(1X,7018.9) 
ID  =  10  ♦  1 

IF(T  .LT.  161  GO  TO  50 
IF(RMAX   .LT,   .001)   CALL  EXIT 


EXC 15400 

EXC15500 

EXC15600 

EXC15700 

EXC15800 

EXC15900 

EXC16000 

EXC  16100 

EXC16200 

EXC 16300 

EXC 16400 

EXC16500 

EXC16600 

EXC16700 

EXC16800 

EXC16900 

EXC17000 

EXC17100 

EXC17200 

EXC17300 

EXC17400 

EXC1750C 

EXC17600 

EXC17700 

EXC17800 

EXC17900 

EXC 18000 

EXC18100 

EXC1820O 

EXC18300 

EXC18400 

EXC18500 

EXC18600 

EXC18700 

EXC18800 

EXC18900 

EXC19000 

EXC19100 

EXC19200 

EXC19300 

EXC19400 

EXC19500 

EXC19600 

EXC19700 

EXC19800 

EXC 19900 

EXC20000 

EXC20100 

EXC20200 

EXC20300 

EXC20400 

EXC20500 

EXC20600 

EXC23700 

EXC20800 


IP1    »    IP1    *    2 

IP2    «    IP2   *    2 

IP3    *    IP3    *   2 

IP*   .    ip*   *    2 

IP5    *    IP5   *    2 

IP6    «    IP6   *    2 

IP7    »    IP7    *    2 

500 

CONTINUE 

CALL    UNDERZ(»ON»> 

CALL    EXIT 

200 

WRITE(6,2000) 

2000 

FORMAT!1    BAD   MATRIX1 » 

CALL    EXIT 

END 

128 


EXC 20900 
EXC21000 
EXC21100 
EXC21200 
EXC21300 
EXC21400 
EXC21500 
EXC21600 
EXC21700 
EXC21800 
EXC21900 
EXC22000 
EXC22100 
EXC22200 


129 
VITA 

Leonard  Andrew  Larsen  was  born  on  April  22,  19UO.  He  attended 
elementary  and  secondary  schools  in  Waukesha,  Wisconsin.  After  graduating 
from  Waukesha  High  School  in  1958,  he  attended  the  University  of  Wisconsin, 
Madison,  where  the  Bachelor  of  Science  degree  in  Physics  and  Mathematics 
Education  was  conferred  in  1962.   Prior  to  attending  the  University  of 
Illinois  in  1967,  he  taught  mathematics  and  physics  at  De  Forest  High  School 
at  De  Forest,  Wisconsin,  for  four  years  and  at  Conant  High  School  in 
Hoffman  Estates,  Illinois,  for  one  year.   He  received  a  Master  of  Science 
degree  in  Mathematics  from  the  University  of  Illinois  in  August  1968.  While 
continuing  his  education  at  the  University  of  Illinois  he  taught  at  Illinois 
State  University  during  the  1969  academic  year.  Beginning  with  the  1971 
academic  year,  he  has  held  the  position  of  assistant  professor  of  computer 
science  at  the  University  of  Wisconsin,  Eau  Claire.  He  is  currently  an 
associate  member  of  the  Association  for  Computing  Machinery  and  Sigma  Xi . 


orm  AEC-427 

(6/68) 
AECM  3201 


U.S.  ATOMIC  ENERGY  COMMISSION 

UNIVERSITY-TYPE  CONTRACTOR'S  RECOMMENDATION   FOR 

DISPOSITION  OF  SCIENTIFIC  AND  TECHNICAL  DOCUMENT 

( See  Instructions  on  Reverse  Side  ) 


AEC  REPORT  NO. 

COO  -ll*69  -0212 


2.    TITLE 


AUTOMATIC  SOLUTIONS  OF  PARTIAL  DIFFERENTIAL  EQUATIONS 


TYPE  OF   DOCUMENT    (Check  one): 

f~1  a.  Scientific  and  technical  report 

[]  b.  Conference  paper  not  to  be  published  in  a  journal: 

Title  of  conference 

Date  of  conference 


Exact  location  of  conference. 
Sponsoring  organization 


|X]  c.  Other  (Specify)       Thesis 


RECOMMENDED  ANNOUNCEMENT  AND  DISTRIBUTION    (Check  one): 

0  a.  AEC's  normal  announcement  and  distribution  procedures  may  be  followed. 

□  b.  Make  available  only  within  AEC  and  to  AEC  contractors  and  other  U.S.  Government  agencies  and  their  contractors. 

1  I  c.  Make  no  announcement  or  distrubution. 


REASON    FOR    RECOMMENDED   RESTRICTIONS: 


SUBMITTED   BY:      NAME    AND   POSITION    (Please  print  or  type) 

C.  W.  Gear,  Professor 

and  Principal  Investigator 


Organization 


Department  of  Computer  Science 
University  of  Illinois 
Urbana,  Illinois  6l801 


Signature 


^KM^d^f^-- 


Date 


October  1972 


FOR   AEC   USE   ONLY 

AEC  CONTRACT  ADMINISTRATOR'S  COMMENTS,   IF    ANY,  ON   ABOVE    ANNOUNCEMENT  AND   DISTRIBUTION 
RECOMMENDATION: 


PATENT  CLEARANCE: 


□  a.  AEC  patent  clearance  has  been  granted  by  responsible  AEC  patent  group. 

□  b.    Report  has  been  sent  to  responsible  AEC  patent  group  for  clearance. 
I    I  c.  Patent  clearance  not  required. 


BIBLIOGRAPHIC  DATA 
SHEET 

4.  Title  and  Subtitle 


1.   Report  No. 


UIUCDCS-R-72-5U6 


AUTOMATIC  SOLUTIONS  OF  PARTIAL  DIFFERENTIAL  EQUATIONS 


7.  Author(s) 


Leonard  Andrew  Lars en 


9.  Performing  Organization  Name  and  Address 

Department  of  Computer  Science 
University  of  Illinois 
Urbana,  Illinois  6l801 

12.  Sponsoring  Organization  Name  and  Address 


US  AEC  Chicago  Operations  Office 
9800  South  Cass  Avenue 
Argonne,  Illinois 


15.  Supplementary  Notes 


3.  Recipient's  Accession  No. 


5.  Report  Date 

October   1972 


8.   Performing  Organization  Rept. 
No. 


10.  Project/Task/Work  Unit  No. 

Ili69-Gear 


11.  Contract /Grant  No. 

US   AEC  AT(ll-l)lU69 


13.  Type  of  Report  &  Period 
Covered 

Thesis    (Ph.D.) 


14. 


16.  Abstracts 


The  problem  of  the  present  paper  is  to  look  for  an  automatic  method 
that  can  be  applied  to  some  subset  of  partial  differential  equations.   The 
primary  value  of  a  program  involving  an  automatic  method  would  be  to 
provide  the  person  who  needs  only  a  few  runs  with  a  particular  type  of 
equation  the  opportunity  to  find  a  solution,  within  acceptable  tolerances, 
without  having  to  write  and  debug  a  specialized  program. 


17.  Key  Words  and  Document  Analysis.     17a.  Descriptors 

partial  differential  equations 
automatic  method 
spatial  discretization 


7b.   Identifiers/Open-Ended  Terms 


7e.  COSAT1  Field/Group 


8.  Availability  Statement 


unlimited  distribution 


ORM   N  TIS-3B    (  10-70) 


19.  Security  Class  (This 
Report) 

VflCLASSIFIED 


20.  Security  Class  (This 
Page 

UNCLASSIFIED 


21-  No.  of  Pages 

129 


22.   Price 


USCOMM-DC    40329-P7  1 


